AD-A070  086 


UNCLASSIFIED 

|of| 

HSrOOB* 


ARMY  MATERIALS  AND  MECHANICS  RESEARCH  CENTER  WATERTOWN  MA  F/6  9/2 
MMCLIBJ  A COMPUTER  CODE  FOR  STRESS  ANALYSIS  USING  THE  MODIFIED  — ETC<U) 
MAR  79  P G TRACY 

AMMR6-TR-79-13  ML 


^MMCLIB:  A COMPUTER  CODE  FOR 


MODIFIED 


l>  MAPPING  -COLLOCATION 


March  1979 


Approved  for  public  release;  distribution  unlimited. 


ARMY  MATERIALS  AND  MECHANICS  RESEARCH  CENTER 
Watertown,  Massachusetts  02172 


The  findings  in  this  report  are  not  to  be  construed  as  an  official 
Department  of  the  Army  position,  unless  to  designated  by  other 
authorized  documents. 

Mention  of  any  trade  names  or  manufacturers  in  this  report 
dtall  not  be  construed  as  advertising  nor  as  an  official 
indorsement  or  approval  of  such  products  or  companies  by 
the  United  States  Government 


' ‘ 1 


DISCLAIMER  NOTICE 


THIS  DOCUMENT  IS  BEST  QUALITY 
PRACTICABLE.  THE  COPY  FURNISHED 
TO  DDC  CONTAINED  A SIGNIFICANT 
NUMBER  OF  PAGES  WHICH  DO  NOT 
REPRODUCE  LEGIBLY. 


'TT7TFT 


UNCLASSIFIED 

||TV  CLASSIFICATION  OF  THIS  PAGE  fit—  Pit  Inf'td) 

J REPORT  DOCUMENTATION  PAGE 


|2.  COVT  ACCESSION  NO. 


/AMMRG-TR-79-13  / j 

s.  iith  7»nSm?«)  — 

^A ^COMPUTE RECODE  FOR^SJRESS^NALYSIS 

. JJSINGTHE  MODIFIED  MAPPING-COLLOCATION  METHOD. 
\ <^  %.  s. 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 
1.  RECIPIENT'S  CATALOG  NUMBER 


rvPE  or  report  a period  covered 


Final  jtfe-pmtt 


SRT  NUMBER 


auTHOAf.j 


“y  Peter  G.^  Tracy  j 

».  PERFORMING  ORGANIZATION  NAME  ANO  ADDRESS 

Army  Materials  and  Mechanics  Research  Center 
Watertown,  Massachusetts  02172  — 

DRXMR-  TM 

• I.  CONTROLLING  OFFICE  NAME  ANO  AO  DREES  i 

U.  S.  Army  Materiel  Development  and  Readines^ 
Command,  Alexandria,  Virginia  22333 


•.  CONTRACT  OR  GRANT  NUMBERED 


io.  program  eleient.  project.  task 

AREA  • WORKBNIT  IINMBERI  -'I 

0/A  Project:  A1T1611^1A91A/ 


Mar*  1§79  / ' 


MONITORING  AGENCY  NAME  4 ADORESV"  Jllloront  Irom  Controlling  Ollico)  I H SECURITY  CLASS,  (ol  thl»  roport) 


Unclassified 

»Sa.  DECLASSIFICATION ’DOWNGRADING 
SCHEDULE 


IS.  DISTRIBUTION  STATEMENT  (ol  tAla  Koporl) 


Approved  for  public  release;  distribution  unlimited. 


I 17.  DISTRIBUTION  STATEMENT  (ol  tho  obatrmct  ontorod  In  Block  20.  II  dlHoronl  Irom  Roport) 


[ i§.  supplementary  NOTES 


I It-  KEY  WOROS  (Contlmto  on  roworoo  aido  II  necssaary  and  Identity  by  block  number) 


Computer  programs 
Stress  analysis 
Complex  variables 


Elasticity 
Numerical  techniques 
Boundary  value  problems 


1 20.  ABSTRACT  (Continue  on  rtmia  aid*  II  neceaaerymnd  I don  Illy  byblock  numbor) 


(SEE  REVERSE  SIDE) 


79  06  19  002 


oo  1473  EDITION  OF  I NOV  SS  IS  OBSOLETE 


4 *o  3 id 5 


UNCLASSIFIED 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  l 


Block  No.  20 


ABSTRACT 


— -^*A  computer  code  for  the  solution  of  a wide  variety  of  plane  isotropic 
stress  problems  using  the  modified  mapping-collocation  method  is  presented 
The  mathematical  approach,  capabilities,  and  structure  of  the  code  are  ex- 
amined. Finally,  a sample  problem  is  solved  to  illustrate  the  ease  with 
which  stress  solutions  can  be  obtained  utilizing  the  code. 


SCCUHITV  CL  AI14FIC  ATION  OF  TNI*  FAOKfMM  Oma  Cnrmtf) 


CONTENTS 


Page 


INTRODUCTION 1 

COMPLEX  VARIABLE  FORMULATION 1 

THE  MODIFIED  MAPPING-COLLOCATION  METHOD  5 

THE  CODE  hWCLIB 7 

SAMPLE  PROBLEM 10 

CONCLUSIONS 11 

APPENDIX  A.  IMPORTANT  MODULES  OF  NMCLIB 15 

APPENDIX  B.  SOURCE  LISTINGS  FOR  MMCLIB  17 


I 


j 


INTRODUCTION 

The  stress  analysis  for  plane  problems  in  elasticity  has  been  handled  using 
a variety  of  methods,  which  include  finite  element  and  finite  difference  methods, 
integral  equations,  integral  transform  techniques,  complex  variable  methods,  and 
boundary  collocation.  One  particularly  versatile  technique  is  the  modified 
mapping-collocation  (MMC)  method1  and  the  MMC  method  with  partitioning2  developed 
at  AMMRC  by  a team  of  investigators  led  by  0.  L.  Bowie.  This  method  combines  the 
complex  variable  and  conformal  mapping  techniques  of  Muskhelishvili3  with  boundary 
collocation  methods.  Finite  element  representations  have  been  introduced  into 
problems  analyzed  by  the  MMC  method  by  Freese,4  resulting  in  greater  potential 
versatility  for  the  MMC  method. 

Previous  exploitation  of  the  hWC  method  to  plane  problems  of  practical  inter- 
est has  been  hampered  by  the  time-consuming  task  of  algebraic  manipulation  and 
subsequent  computer  programming  of  the  equations  used  in  a given  problem.  Since 
the  equations  used  are  similar  in  each  problem,  the  computer  code  MMCLIB  (Modified 
Mapping-Collocation  LIBrary)  has  been  written  to  eliminate  the  algebraic  manipula- 
tion effort  and  greatly  minimize  the  computer  programming  task.  MMCLIB  consists 
of  a library  of  ninety  Fortran  subroutines  and  functions  which,  when  combined  in 
a Fortran  main  program  written  by  the  user,  have  the  ability  to  quickly  provide 
stress  solutions  for  plane  problems  while  retaining  the  versatility  inherent  in 
the  MMC  method.  Furthermore,  addition  of  new  mapping  functions  and  series  rep- 
resentations can  be  accomplished  with  a minimum  of  effort,  which  allows  the  ex- 
pansion of  the  code's  capabilities.  The  versatility  of  the  MMC  method  has  been 
demonstrated  in  the  solution  of  a wide  variety  of  stress  problems  such  as  sheets 
with  internal  flaws,  i.e.,  cracks  and  ellipses,  sheets  with  edge  cracks  and  ellip- 
tical notches,  sheets  with  cracks  emanating  from  elliptical  internal  flaws  and 
edge  notches,  and  circular  rings  with  edge  cracks. 

The  current  report  is  not  intended  to  be  a user's  manual  but  is  intended  to 
describe  the  capabilities  o¥~ the  code  MMCLIB  in  the  solution  of  plane  isotropic 
stress  problems.  Information  on  the  use  of  MMCLIB  can  be  obtained  from  the 
author . 


COMPLEX  VARIABLE  FORMULATION 


The  complex  variable  approach  provides  the  mathematical  formulation  on  which 
the  MMC  technique  is  based  and  results  relating  to  the  MMC  method  will  be  discussed 
in  this  section.  For  a complete  treatment  of  the  subject  the  reader  is  referred  to 
Muskhelishvili. 3 A well-known  method  of  describing  the  stresses  in  an  isotropic 


BOWIE,  O.  L.,  and  NEAL,  D.  M.  A Modified  Mapping  Collocation  Technique  for  Accurate  Calculation  of  Stress  Intensity  Factors. 
Int.  J.  Fracture  Mech.,  v.  6,  1970,  p.  199-206. 

BOWIE,  O.  L.,  FREESE,  C.  E.,  and  NEAL,  D.  M.  Solutions  of  Plane  Problems  of  Elasticity  Utilizing  Partitioning  Concepts. 

J.  Appl.  Mech.,  v.  95,  1970,  p.  767-772. 

MUSKHELISHVILI,  N.  I.  Some  Basic  Problems  of  the  Mathematical  Theory  of  Elasticity.  Noordhoff,  Groningen,  Holland,  1953. 
FREESE,  C.  E.  Collocation  and  Finite  Elements  ■ A Combined  Method.  Army  Materials  and  Mechanics  Research  Center,  AMMRC 
TR  73-28,  1973. 


elastic  two-dimensional  body  is  through  the  use  of  the  Airy  stress  function,  U(x,y). 
Specifically,  equilibrium  and  compatibility  conditions  will  be  satisfied  if 


f 


v2v2  U(x,y)  * 0 (l) 

with  the  stresses  defined  by 
ox  - [32U(x,y)]/3y2 
oy  - [32U(x,y)]/3x2 

xxy  - [-32U(x,y)]/3x3y  (2) 

where  V2  is  the  two-dimensional  Laplacian  operator.  In  the  Muskhelishvili  formula- 
tion, the  coordinates  of  the  problem  are  described  using  the  complex  variable  z ■ 
x ♦ iy  and  the  z-plane  is  referred  to  as  the  physical  plane.  The  Airy  stress  func- 
tion can  then  be  represented  using  two  analytic  functions  $z(z)  and  ^z(z)  by 

U(x,y)  * Re  [T  $z(z)  + X2(z)]  C3) 

where  i|>z(z)  = d Xz/dz,  Re  means  the  real  part,  and  bars  complex  conjugation.  The 
stresses  and  displacements  in  terms  of  $z(z)  and  4>z(z)  are 

CTX  + <Ty  = 4 Re  [$z(z)  ] 

ay  - ax  + 2i  txy  = 2 [T  (z)  + rz( z)  ] (4) 

2 6 (u  ♦ iv)  » k<|>z(z)  - zJiJzT  - 0zCz)  (5) 

with  G = E/2  (1+v),  and  k = (3-v)/(l+v)  (plane  stress),  or  k = 3-4v  (plane  strain), 
E being  Young's  modulus  and  v being  Poisson's  ratio. 

An  additional  physical  quantity  which  is  used  frequently  in  applying  the  MMC 
technique  is  the  force  acting  on  an  arbitrary  arc  of  material.  If  Xnds  and  Ynds 
are  x-  and  y-components  of  the  force  acting  on  an  element  of  the  arc  ds  (see 
Figure  1),  it  can  be  shown  that 

(X„  ♦ i Yn)ds  = -i  d (*z(z)  ♦ zjijz)  ♦ r^Z)} 


Figure  1.  Forces  acting  on  an  arc  of  material 


The  resultant  moment  about  the  origin  of  coordinates  can  be  expressed  as 

M = f{x  Yn  - y Xn)  ds 

- Re  [Xz(z)  - z i|»z(z)  - zz  ♦•(z)]J  (7) 

where  A and  B are  the  endpoints  of  the  arc. 

The  formulation  of  the  two-dimensional  theory  of  elasticity  in  terms  of  ana- 
lytic functions  permits  the  use  of  conformal  mapping  techniques,  providing  the 
means  to  introduce  curvilinear  coordinates  to  describe  the  geometry  in  the  physi- 
cal plane.  Use  of  conformal  mappings  leads  to  the  introduction  of  the  c or  param- 
eter plane  with  the  conformal  transformation 

z = u(C).  (8) 

The  functions 


4>CO  = *z  (u>(0) 
<KO  = ¥z  (u(C)) 


(9) 


are  therefore  analytic  functions  of  e since  <frz(z)  and  <|>z(z)  are  analytic  functions 
of  z and  w(0  = z is  a conformal  mapping.  Substituting  <|>(0  and  <|>(0  for  <t>z(z) 
and  ipziz)  and  using  relationships  such  as  <fr£(z)  = ' (^)/aj'  (c)  > Equations  4,  5,  6, 

and  7 become 


ox  + Oy  = 4 Re 


IV  CO 
pTcTj 


- o + 2i  t =2  (*••(£)  ,-  + £110 

°x  + Txy  * (oj'  CO  i2  ^ UJ  «'  CO  J »'  CO 


2G  (u+iv)  = k $C0  - 


w(0 


u'CO 


*’C0  - <ko 


+ if2  = 4>(0 


. toCO 


<*>’  CO 


4-'  CO  + £(0 


m = Re  [XCO  - <o(0<KO  - WCO  «(0  ♦•CO/w'COlX- 


CIO) 

C11) 

(12) 

(13) 


An  extremely  useful  result  of  the  work  of  Muskhelishvili  is  the  application 
of  the  analytic  continuation  arguments  summarized  below.  Let  S+  denote  the  region 
above  the  real  axis  (y  > 0)  and  S“  the  region  below  the  real  axis  (y  < 0)  and  fur- 
ther assume  that  the  elastic  body  occupies  the  region  S+  (see  Figure  2) . The 
functions  $z(z)  and  ij>z(z)  are  thus  defined  in  S+.  The  functions  4>z(z)  shall  be 
continued  into  S"  by  the  following  definition: 

♦ z(z)  * -zf^Cz)  - iTz(z)  for  zeS",  (14) 
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Figure  2.  Analytic  continuation. 


with  the  notation  F(z)  = f(l) . Noting  that  if  zeS 
plex  conjugate  of  Equation  14  one  finds 


then  zeS  and  taking  the  com 


Substitution  of  Equation  15  into  Equation  6 gives 
fx  + if2  = <t»z(z)  - <|>zCz)  ♦ (z-z)  *£(z)  . 


Evaluation  of  this  equation  for  z with  Im(z)  * 0,  where  Im  means  imaginary  part, 
gives  fi  + if 2 = 0.  Thus,  the  real  axis  has  been  defined  as  load-free,  while 
i|)z(z)  has  been  eliminated  from  the  analysis  leaving  only  the  analytic  function 
4>z(z)  to  be  determined.  If  Equation  15  is  also  substituted  into  Equations  4 and 
5 the  following  results: 


2 G(u+iv)  = k*z(z)  ♦ *z(z)  ♦ (T-z)  ♦•(z) 


Similarly  the  function  $_(z)  can  be  continued  from  the  exterior  of  the  unit 
circle  (S+)  to  the  interior  of  the  unit  circle  (S')  using  the  following  definition 
due  to  Kartzivadze5 


This  leads  to  expressions  similar  to  Equations  17  and  18 


5.  KARTZIVADZE,  I.  N.  The  Fundamental  Problems  of  the  Theory  of  Elasticity  for  the  Elastic  Circle.  Comp.  Rend,  de  L'acad.  Sc. 
de  1TJRSS,  v.  20,  1943,  p.  95-104. 


(21) 


r 


2 G(u*iv)  = k$z(z)  ♦ ♦2Cl/z)  ♦ (1/z  - z)  ♦i(z) 
fi  ♦ if2  = <>z(z)  - ♦jfl/T)  ♦ (z  - 1/z)  ♦i(z)  . 


The  evaluation  of  Equation  21  on  the  unit  circle  shows  that  traction- free 
conditions  must  exist  on  the  unit  circle.  It  should  be  noted  that  the  definitions 
in  Equations  19  and  13  for  the  analytic  continuation  of  $z(z)  from  S+  to  S"  will 
need  slight  modification  when  conformal  mapping  is  included  in  the  analysis  and 
♦ (C)  = ♦z(io(C))  is  used  due  to  the  presence  of  the  mapping  function  in  the  force 
Equation  12. 


THE  MODIFIED  MAPPING-COLLOCATION  METHOD 

The  modified  mapping -col location  0#C)  method  is  a technique  for  two- 
dimensional  elastic  stress  analysis  which  combines  the  complex  variable  formula- 
tion and  conformal  mapping  arguments  of  Muskhelishvili3  and  the  boundary  colloca- 
tion technique,  retaining  the  attractive  features  of  each.  This  method  was  first 
proposed  by  Bowie  and  Seal,1  and  was  refined  with  the  introduction  of  a partition- 
ing plan  by  Bowie,  Freese,  and  Neal.2  The  >WC  method  has  subsequently  been  applied 
to  a wide  variety  of  problems,6-11  illustrating  its  versatility.  Refer  to  Figure  3 
for  illustrations  of  some  of  the  problems  solved. 

In  the  MMC  technique  a simple  mapping  function  as  in  Equation  7 is  chosen  to 
describe  critical  portions  of  the  physical  geometry,  e.g.,  cracks,  notches,  cutouts, 
etc.  Remaining  portions  of  the  geometry  can  be  described  by  inversion  of  the  plane 
mapping  function  from  the  physical  to  the  parameter  plane.  Thus,  the  forms  of  the 
mapping  functions  are  algebraically  simple  since  the  entire  geometry  need  not  be 
described  by  a single  function,  and  greater  versatility  is  obtained  since  a unique 
mapping  function  need  not  be  obtained  for  each  problem  geometry.  The  specifica- 
tion in  the  parameter  plane  of  those  portions  of  the  boundary  not  described  by  the 
mapping  function  may  be  algebraically  complex,  making  direct  specification  of 
boundary  conditions  extremely  difficult.  For  this  reason  a boundary  collocation 
scheme  is  introduced. 

In  order  to  obtain  $(£)  and  i/»(C)  for  the  problem  under  consideration  a rep- 
resentation for  each  function  must  be  chosen.  Generally,  when  analyticity  is 
assured,  power  and  Laurent  series  are  chosen  for  representation,  and  the  unknown 
coefficients  of  the  series  are  determined  using  the  boundary  collocation  arguments. 
The  forms  for  <HO  and  <Ji(0  are  substituted  into  Equations  10  to  13,  and  equations 


6.  BOWIE,  O.  L.  Methods  of  Analysis  and  Solutions  of  Crack  Problems  in  Mechanics  of  Fracture,  v.  1,  G.  C.  Sih,  cd.,  Noordhoff, 
Leden,  1973. 

7.  BOWIE,  O.  L.,  and  FREESE,  C.  E.  Analysis  of  Notches  Using  Conformal  Mapping  in  Mechanics  of  Fracture,  v.  S,  G.  C.  Sih,  cd., 
Noordhoff,  Leden,  1978. 

8.  BOWIE,  O.  L.,  and  FREESE,  C.  E.  On  the  'Overlapping'  Problem  in  Crack  Analysis.  Engineering  Fracture  Mechanics,  v.  8,  1976, 
p.  373-379. 

9.  BOWIE,  O.  L.,  and  FREESE,  C.  E.  Plastic  Analysis  for  a Radial  Crack  in  a Circular  Ring.  Army  Materials  and  Mechanics 
Research  Center,  AMMRC  MS  70-3,  1970. 

10.  TRACY,  P.  G.  Plastic  Analysis  of  Radial  Cracks  Pmanating  from  the  Outer  and  Inner  Surfaces  of  a Circular  Ring.  Engineering 
Fracture  Mechanics,  v.  II,  1979,  p.  291. 

1 1.  TRACY,  P.  G.  Analysis  of  a Radial  Crack  in  a Circular  Ring  Segment.  Engineering  Fracture  Mechanics,  v.  7,  1975,  p.  253-260. 
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can  be  written  in  terms  of  the  unknown  coefficients  of  <KO  and  In  situa- 

tions where  the  representations  for  $(c)  and  i/»(c)  are  complete,  collocation  of  fj 
and  f2  (Equation  12)  is  usually  sufficient.  In  a region  with  a crack  tip  near  a 
boundary,  error  in  matching  fj  and  f2  may  result,  indicating  the  need  for  collocat- 
ing stresses.  Equation  10,  as  well  as  fj  and  f2,  particularly  on  the  portion  of 
the  boundary  near  the  crack  tip.  In  cases  where  the  representation  may  not  be 
complete,  the  collocation  of  M,  Equation  13,  as  well  as  forces  and  stresses  may 
be  required.  A set  of  points  on  the  boundaries  is  chosen  at  which  to  write  equa- 
tions such  that  the  number  of  equations  generated  is  greater  than  the  number  of 
unknown  coefficients  by  at  least  a factor  of  two.  A least-squares  minimization 
of  error  is  used  to  generate  a linear  system  of  equations  in  the  unknown  coeffi- 
cients; the  coefficients  can  then  be  determined  by  solving  the  system  of  equations. 

The  partitioning  concepts  introduced  by  Bowie,  Freese,  and  Neal2  extend  the 
flexibility  of  the  MMC  technique  to  allow  solutions  for  problems  which  are  awk- 
ward with  a single  series  representation  for  the  stress  functions.  Such  problems 
occur  when,  for  example,  singularities  occur  near  the  circles  of  convergence  for 
the  stress  functions  due  to  areas  of  high  stress  concentrations.  To  implement 
the  partitioning  plan,  the  geometry  is  subdivided  into  two  or  more  regions.  The 
functions  $2(z)  and  ^z(z)  are  defined  in  a subregion  i by  <J>z,i(2)  and  i|/2>i(z), 
along  with  a mapping  for  each  region,  if  appropriate.  Boundary  conditions  (Equa 
tions  10  to  13)  are  enforced  along  the  boundary  using  the  appropriate  $2  ^(z)  and 
<|>z  i(z)  for  the  collocation  point  at  which  equations  are  written.  Additionally, 
at’points  on  the  common  boundaries  of  zones  u+iv  and  fx  + if2  are  equated  to  as- 
sure continuity  of  <j>2(z)  and  i|/z(z).  Then  the  least-squares  procedure  is  invoked 
and  the  resulting  equations  solved  to  obtain  the  function  <|>z  j(z)  and  \\>z  j(z)  in 
each  subregion. 
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Once  the  solution  is  obtained  it  is  necessary  to  assess  the  degree  to  which 
the  boundary  conditions  are  satisfied  by  $z(z)  and  *z(z)  as  obtained.  This  is 
done  in  two  ways.  First,  the  solution  is  recomputed,  varying  the  number  and  lo- 
cation of  points  at  which  equations  are  written  and  the  number  of  unknown  coeffi- 
cients to  confirm  that  a wide  range  of  these  choices  will  yield  the  same  results, 
thus  assuring  a convergent  solution.  As  a second  check,  the  boundary  conditions 
are  evaluated  using  the  computed  $z(z)  and  tfiz( z)  and  compared  to  the  values  re- 
quired for  solution.  Good  agreement  will  assure  correct  results. 

THE  CODE  MMCLIB 

The  computer  code  MMCLIB  has  been  written  to  enhance  implementation  of  the 
KWC  method.  To  retain  maximum  flexibility  in  the  variety  of  problems  to  be  solved, 
the  code  is  constructed  in  a modular  fashion.  The  user  is  required  to  write  a 
Fortran  main  program  which  usually  is  30  to  100  lines  in  length  and  uses  the  mod- 
ules (Fortran  subroutines  and  functions)  of  MMCLIB.  In  addition,  the  modular  con- 
struction allows  the  implementation  of  additional  mapping  functions  and  series 
representation  types  as  one  need  program  only  the  function  or  representation  and 
its  first  two  derivatives  into  the  code.  Let  us  now  consider  the  solution  process, 
which  is  illustrated  in  Figure  4,  and  examine  the  features  of  MCLIB. 


Figure  4.  Block  diagram  of  solution  process  using  MMC  method. 
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The  problem  geometry  under  consideration  must  first  be  subdivided  into  one  or 
more  subregions  or  zones  with  reflections  as  in  Equation  14  or  Equation  19,  if 
appropriate.  This  will  define  the  number  and  type  (either  $z>i(z)  or  i|>z  i(z))  of 
functions  needed  to  complete  the  solution.  Each  such  function  is  defined  by  zone 
in  MMCLIB  by  a subroutine  call.  At  present  two  functional  representations  are 
available: 

fi(W)  = £ - «o)n  (22) 

n=-N 


and 


f2(W) 


an(Wn  - W0)n 
u W+i 


(23) 


where  an  is  an  unknown  complex  coefficient,  W is  the  (complex)  parametric- 
coordinate  or  z- coordinate  in  the  case  of  no  mapping,  M and  N are  nonnegative 
integers,  and  WQ  is  the  expansion  point.  A function  is  defined  by  specifying  WQ, 
M,  N,  the  series  representation,  the  material  parameters  of  that  zone,  and  any 
real  or  imaginary  parts  of  an  for  even  or  odd  powers  of  n Which  are  assumed  to 
be  zero  due  to  the  stress  symmetries  of  the  problem.  For  each  zone  a mapping 
function  may  be  chosen  from  the  following  two  presently  available 

coi  (C)  = [(ai+bjK/2  + (aj-bjK'Vzie*6!  + C0,l  (24) 

and 

u>2 (C)  = i[b2cosh(ic)  - a2sinh(i?) ]ei02  + Co, 2-  (2S) 


The  action  of  these  two  mapping  functions  is  shown  in  Figures  5 and  6.  In  addi- 
tion, MMCLIB  allows  mapping  functions  which  are  compositions  of  <di(c)  and  w2(C), 
i.e.. 


<*>(?)  = b»i[<»j(0] 


(26) 


where  i and  j vary  from  1 to  the  number  of  mappings  currently  available. 


One  of  the  most  tedious  aspects  of  problem  solving  using  the  WC  method  has 
been  the  need  to  algebraically  reduce  equations  of  the  type  Equations  10  to  13, 
17,  18,  20,  and  21  with  the  substitutions  for  <|>z(z)  and  i|iz(z)  or  *(C)»  *(C)  of 
representations  of  the  form  of  Equations  22  and  23  to  the  form 


[“nr  frJn)P0  * <ni  fr{#)00]  = R® 


u+iv 
fl  ♦ if2 

°y  “ °x  + 2:*-  Txy 
4 ♦£(*) 
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and 


t^nr  fi£nJ(W)  ♦ ani  f^W]  = In. 


*1  ♦ if2 


0X  + X Xy 


where  an  - anr  + i a^,  W is  either  ; or  z and  the  f's  are  functions  of  the  posi- 
tion of  the  point  at  which  equations  are  being  written.  Collocation  of  forms  for 
functions  such  as  a<j>z(z),  a<j>z(z),  (z),  a<t>,(T),  <*4>z(z) , etc.,  are  each  contained 

in  a subroutine,  and  equations  are  generated  by  successive  calls  to  these  subrou- 
tines. Note  that  a is  some  complex  number.  The  form  assumed  for  each  function  is 


£ an  fnW  (28) 

n 

where  presently  fn(W)  is  either  of  Equation  22  or  23.  Other  forms  for  fn(W)  can 
be  added  by  simply  inserting  programming  for  fn(W),  [dfn(W)/dW],  and  [d2fn(W)/d2W] 
into  the  three  appropriate  subroutines.  Many  of  the  forms  for  the  mapping  func- 
tions such  as  oj(c),  w(C),  <*>'({),  &)''({),  w({),  etc.,  are  also  contained  in  modules 
for  use  in  generation  of  equations.  It  should  be  noted  that  if  &>i({),  <4(c),  a|'  (c) 
are  programmed  into  the  three  appropriate  subroutines,  the  mapping  function 
can  be  included  for  analysis.  Certain  commonly  used  equations  for  problem  solving 
have  been  programmed  into  modules  using  previously  described  subroutines.  Thus, 
equations  for  forces,  displacements,  and  stresses  for  real  axis  reflections,  unit 
circle  reflections,  and  no  reflection  are  contained  in  modules.  MMCLIB,  therefore, 

I C -plane 


Figure  5.  The  mapping  function  z = (f ). 


I {-plane 


MclK 


L <*_ 


Figure  6.  “he  mapping  function 
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eliminates  most  of  the  repetitious  equation  derivation,  allowing  the  user  to  con- 
centrate on  modelling  the  problem  at  hand  without  spending  time  in  algebraic  manip- 
ulation of  equations. 

After  points  are  chosen  for  collocation  and  the  equations  in  the  unknown  co- 
efficients are  generated  and  solved,  the  results  of  the  computations  can  be  exam- 
ined. Any  physical  quantity  can  then  be  calculated  using  the  computed  coefficients 
and  the  appropriate  modules  to  generate  equations  for  the  quantities  of  interest. 
The  more  common  physical  quantities  of  stress,  displacements,  and  forces  have  been 
programmed  into  modules  to  facilitate  the  acquisition  of  results.  Subroutines  are 
also  available  to  print  out  the  computed  coefficients  and  the  input  information 
about  the  various  series  in  use.  See  Appendix  A for  the  names  and  functions  of 
the  important  modules  of  MMCLIB  and  Appendix  B for  the  Fortran  source  listings  for 
MMCLIB. 

SAMPLE  PROBLEM 

To  illustrate  the  ease  with  which  solutions  can  be  obtained  using  NMCLIB,  the 
problem  of  an  elliptical  edge  notch  in  a semi-infinite  sheet  under  tension  was 
solved  (Figure  7).  This  problem  was  studied  by  Bowie  and  Freese7  and  their  method 
of  solution  was  duplicated  using  FWCLIB.  The  mapping  function  used  for  zone  I is 

z = (oU)  = i [(a+b)e"i?/2  - (a-b)e-iC/2]  (29) 

which  is  a composition  of  mapping  functions  of  the  type  in  Equations  24  and  25. 


Figure  7.  Sample  problem:  Elliptical  edge  notch  in 
semi-infinite  sheet  under  remote  tension. 


A a B 


<*>(c)  = (io2  CO ) 


ai  = a bi  = b 


in  Equation  24  and 
a2  = b2  = 1 


el  = 0 Co  l = 0 


0 = 0 Co, 2 * 0 


in  Equation  25.  This  mapping  takes  the  real  axis  in  the  c -plane  such  that  - ir/2  < 
Re  c 5 it/ 2 into  the  elliptical  arc  AB  in  the  z-plane.  In  order  to  ensure  that  AB 
is  traction- free,  a reflection  across  the  real  axis  in  the  c-plane  is  used  with  the 
following  representation  for  $z(z)  = ♦z(u(C))  * $i(c)  in  zone  I 


(31) 


Ni 

<*>l(C)  * 5Z  ta2n-l  52n_1  + i a2n  C2n] 
n=l 

where  the  a's  are  real,  ensuring  stress  symmetry  about  the  physical  imaginary  axis. 
In  zone  II  there  is  no  mapping  function  and  a reflection  is  used  to  ensure  traction- 
free  conditions  on  the  real  azis.  The  representation  for  <t>z(z)  is 

N2 

4>Il(z)  = (1/4)  + 5Z  ta-2n- 1 z'2n_1  + i a-2n  z_2n]  (32) 

n=0 

ensuring  ox  = aQ  for  large  z and  the  proper  stress  symmetry  about  the  imaginary 
axis.  Along  the  segment  BD  the  condition  fj  + if2  = 0 is  collocated,  and  the  con- 
dition of  continuity  of  f^  + if2  and  u+iv  are  collocated  on  the  circular  arc  CD. 

The  Fortran  main  program  used  in  the  solution  of  this  problem  is  shown  in 
Figure  8 to  illustrate  the  magnitude  of  the  computer  programming  effort  necessary 
to  carry  out  a solution  using  the  MMC  method  with  the  code  MMCLIB.  Lines  8 to  12 
define  the  series  and  mapping  functions  to  be  used  and  lines  16  to  20  and  22  to 
30  find  points  and  collocate  equations  on  the  curves  BD  and  DE,  respectively.  The 
printing  of  results  of  the  computations  are  programmed  into  lines  34  to  40.  It 
can  be  seen  that  programming  effort  to  solve  a problem  using  the  MMC  method  is 
minimal  with  WCLIB. 

The  computer  generated  output  for  this  problem  is  shown  in  Figure  9.  The 
maximum  stress  as  computed  by  MMCLIB  compares  well  with  that  computed  by  Bowie 
and  Freese.  The  error  calculation  in  this  example  is  typical  of  many  problems. 

Note  that  the  right-hand  sides  of  the  equations  generated  are  of  the  order  of 
magnitude  of  the  coordinates  (i.e.,  unity).  The  amount  of  error  compared  to  the 
order  of  magnitude  of  the  values  desired  as  shown  in  the  first  two  columns  of  the 
error  computation  are  usually  considered  adequate  for  accurate  results.  The  third 
and  fourth  columns  and  fifth  and  sixth  columns  are  the  coordinates  in  zones  one 
and  two,  respectively,  at  which  errors  were  checked. 


CONCLUSIONS 

The  computer  code  MMCLIB  should  permit  a savings  in  the  effort  involved  in 
obtaining  stress  solutions  to  plane  problems  of  elasticity  using  the  MMC  method 
since  the  effort  will  be  concentrated  on  the  mathematical  and  physical  concepts 
involved  in  solving  the  problems  and  not  on  mechanical  activities  such  as  alge- 
braic manipulation  of  equations.  Thus,  problems  whose  solutions  are  not  concep- 
tually difficult  to  the  experienced  user  of  the  MMC  method  but  which  were  tedious 
to  implement  can  be  solved  almost  routinely  using  MMCLIB.  Furthermore,  research 
into  further  development  of  the  MMC  method  can  be  performed  much  more  easily 
since  most  of  the  analyst's  attention  can  be  focused  on  the  new  features  being 
developed. 
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DIMENSION  X(50) ,Y(50) ,XP(50) ,YP(50),F1(50),F2(50),U(50) ,V(50) . 

* GPAR  ( 5 ) 

COMMON  /MAT/  G(5),ETAM(5) 

COMMON  /P I FAC/  PI .HLFPl .TWOPl 

COMMON  /SOLN/  A(15401),IDIAG(176).PHI(176) 

DATA  GPAR  /5»0 . / 

CAtL  INIT 

RE AO  1000, IPOS, INEG.NPT.NST.AA.BB 

CALL  SER I NF ( 1 , 1 , 1 , 0 , IPOS.O, 1001 , 1 ,0. .0. , 1 .E+7, .3) 

CALL  SERINF(2. 1 , 1 .INEG.O. 1 , 1001 . 1 ,0. . 0. , 1 . E+7, .3) 

CALL  MAPINF( 1,2,1.,1.,0.,0.,0..1,AA,BB,0..0.,0.) 

RZ«AMAX1 (AA,BB)*2. 

I ORR-O 

100  CALL  P0INTF(1  . AA.O. ,RZ.O.  .NPT.1 , 1 , 1 .GPAR.1 , 2 , X , Y , XP, YP) 

DO  110  1=1 , NPT 
FI ( I ) =0 . 

110  F2( I ) =0 . 

CALL  FDCOHXP.YP.NPT  ,1  , 1 ,2,F1  .F2.IERR) 

GPAR ( 1 ) =RZ 

CALL  P0INTF(1 ,RZ,0. ,0. , RZ . NST ,0 , 0 . 2 .GPAR . 0 . 2 , X , Y , XP , YP) 

DO  120  1=1 . NST 
F 1 ( I ) =0 . 

F2( I )*Y( I ) 

U(I)=X(I)*(ETAM(1)+1.)*.2S 
120  V( I )=Y( I )*(ETAM( 1 )-3 . )*.25 

CALL  STITCH (1 ,2, 2, 2. 3, 2,0. NST ,X,Y,F1,F2tU,V, I ERR) 

I F ( I ERR . EQ • 1 ) GO  TO  130 

CALL  FNDMAT 

PRINT  1002.AA.BB 
CALL  PRTCOFf 1 2 ) 

CALL  ZERO ( 1 ,2. 1 ,2.3,0) 

CALL  PHIP( .0001  ,0. .4./AA.0.  , 1 ) 

CALL  PLACE ( 1 , 1 , 1 . 1 . 1 . 1 . ) 

CALL  EVAL ( 1 , 1 , 1 , STRMAX , SDUM) 

PRINT  1003, STRMAX 

PRINT  1001 
I ERR  = 1 
NST.6 
NPTalO 
GO  TO  100 

130  PRINT  1004 
STOP 

1000  FORMAT (4I5.2E10.4) 

1001  FORMAT ( 7HOERRORS//4X , 7HERR0R  1.3X.7HERR0R  2.3X.6HC0R  1R.4X, 

* 6HC0R  1I.4X.6HC0R  2R.4X.6HC0R  21) 

1002  F0RMAT(45H1ELLIPTICAL  EDGE  NOTCH  IN  SEMI-INFINITE  SHEET// 

* 4H  A = , F8.4. 10X.3HB  «.F8.4/////) 

1003  FORMAT! ////I 7H  MAXIMUM  STRESS  «.F8.4////) 

1004  FORMAT ( 1 HI ) 

END 


Figure  8.  Main  program  for  sample  problem. 


ELLIPTICAL  EOGE  NOTCH  IN  SEMI-INFINITE  SHEET 
A - .1250  B * .5000 


COEFFICIENTS  IN  PHI.PSI  FUNCTIONS 


ZONE  NO.  2 


PHI  FUNCTION 
SERIES  TYPE  1 

ZONE  NO.  1 


POW 

REAL  COEF 

IMAG  COEF 

PHI 

FUNCTION 

-20 

.00000 

.34927-05 

SERIES  TYPE  1 

-19 

-.19418-04 

.00000 

-18 

.00000 

— . 47040-04 

POW 

REAL  COEF 

IMAG  COEF 

-17 

.63436-04 

.00000 

0 

.00000 

.00000 

-16 

.00000 

.32696-04 

1 

.30070+00 

.00000 

-15 

.17873-04 

.00000 

2 

.00000 

.39187-01 

-14 

.00000 

.47948-04 

3 

-.40076-01 

.00000 

-13 

-.11984-04 

.00000 

4 

.00000 

-.85981-02 

-12 

.00000 

.77529-04 

5 

.38902-02 

.00000 

-11 

-.14141-03 

.00000 

6 

.00000 

-. 16099-03 

-10 

.00000 

-.14408-03 

7 

.22604-03 

.00000 

-9 

.19836-03 

.00000 

8 

.00000 

-.  17401-04 

-8 

.00000 

.28798-03 

9 

.51200-04 

.00000 

-7 

-.67002-03 

.00000 

10 

.00000 

. 1 1 194-04 

"6 

.00000 

-.12855-02 

1 1 

.41939-05 

.00000 

-5 

.34014-02 

.00000 

12 

.00000 

.13837-05 

-4 

.00000 

.63706-02 

13 

.13354-05 

.00000 

-3 

-.19156-01 

.00000 

14 

.00000 

.13379-05 

-2 

.00000 

- . 33425-01 

15 

-.37054-06 

. 00000 

-1 

.16262+00 

.00000 

16 

.00000 

-.36768-07 

0 

.00000 

.22439+00 

17 

.35488-07 

.00000 

18 

.00000 

.55372-07 

19 

-.23934-07 

.00000 

20 

.00000 

-.44828-08 

MAXIMUM  STRESS  - 9.6224 


ERRORS 


ERROR  1 

ERROR  2 

COR  1 R 

COR  11 

COR  2R 

COR  21 

-.0003 

-.0006 

.1571+01 

.1888+00 

-.0003 

-.0000 

.1571+01 

.3640+00 

.0005 

-.0000 

.1571+01 

.5239+00 

.0001 

.0000 

.1571+01 

.6688+00 

-.0003 

.0003 

.1571+01 

.7998+00 

-.0001 

.0001 

.1571+01 

.9187+00 

.000 1 

-.0003 

.1571+01 

.1027+01 

.0000 

-.0001 

.1571+01 

.1126+01 

.0001 

.0004 

.1571+01 

.1217+01 

.1000+01 

.0000 

.0000 

-.0000 

.1571+01 

.1217+01 

.1000+01 

.0000 

.0000 

.0001 

.1287+01 

.1209+01 

.9511+00 

.3090+00 

-.0000 

-.0001 

.1287+01 

.1209+01 

.9511+00 

.3090+00 

-.0001 

-.0001 

.9949+00 

.1185+01 

.8090+00 

.5878+00 

.0000 

.0000 

.9949+00 

.1185+01 

.8090+00 

.5878+00 

.0001 

.0001 

.6866+00 

.1150+01 

.5878+00 

.8090+00 

-.0000 

-.0001 

.6866+00 

.1150+01 

.5878+00 

.8090+00 

-.0000 

-.0001 

.3542+00 

.1114+01 

.3090+00 

.9511+00 

.0000 

.0001 

.3542+00 

.1114+01 

.3090+00 

.9511+00 

-.0000 

.0001 

.2980-07 

.1099+01 

.1589-07 

.1000+01 

.0000 

-.0000 

.2980-07 

.1099+01 

.1589-07 

.1000+01 

Figure  9.  Computer  genereted  output  for  temple  problem. 


APPENDIX  A.  IMPORTANT  MODULES  OF  MMCLIB 


In  this  section  the  more  commonly  used  subroutines  of  NWCLIB  will  be  listed 
and  their  uses  explained.  The  subroutines  are  broken  up  by  function  into  five 
groupings:  Zone  Definition,  Mapping  Functions,  Function  Collocation,  Specialized 
Equation  Generation,  and  Results  Computation. 

The  Zone  Definition  routines  allow  the  definition  of  functions  (either  tz(z) 
or  *z(z))  in  each  zone  to  be  used  in  the  analysis  and  definition  of  mapping  func- 
tions in  each  zone.  In  addition,  the  functional  representations,  material  prop- 
erties, and  any  unknowns  needed  but  not  included  in  the  functional  forms  chosen 
are  defined.  Note  that  all  zone  definition  subroutine  calls  are  preceded  by  a 
call  to  INIT.  Table  A-l  lists  these  routines  and  their  uses. 

The  Mapping  Function  routines  allow  the  calculation  of  coordinates  in  either 
the  parametric  or  physical  plane.  These  routines  are  used  in  determination  of 
coordinates  of  points  in  one  plane  when  the  corresponding  coordinates  in  the 
other  plane  are  known.  The  Mapping  Function  routines  and  the  forms  of  the  map- 
ping function  which  they  compute  are  listed  in  Table  A-2.  Note  that  some  subrou- 
tines return  a number  of  functional  forms  depending  on  the  number  of  mapping 
functions  defined  in  a given  zone. 

The  Function  Collocation  subroutines  are  used  to  write  equations  for  the  un- 
known parameters  of  the  functions.  Several  Function  Collocation  routines  are 
listed  in  Table  A- 3. 

The  Specialized  Equation  Generation  routines  are  used  to  generate  multiple 
equations  for  the  most  commonly  used  physical  quantities.  These  routines  are 
frequently  used  since  in  most  solutions  using  the  MMC  method,  the  equations  col- 
located on  a given  part  of  the  boundary  of  the  configuration  being  studied  are 
of  the  same  type  and  thus  it  is  convenient  to  generate  all  these  equations  with 
one  subroutine  call.  Specialized  Equation  Generation  routines  make  use  of  Func- 
tion Collocation  and  Mapping  Function  routines  to  generate  the  equations  for 
physical  quantities.  The  Specialized  Equation  Generation  routines  are  listed  in 
Table  A-4. 

The  Result  Computation  routines  are  used  to  print  the  coefficients  computed 
in  the  solution  and  physical  quantities  which  can  be  obtained  using  the  solution 
for  the  unknown  coefficients.  The  routines  are  also  used  to  print  the  various 
parameters  input  in  the  zone  information  routines.  The  Result  Computation  sub- 
routines are  listed  in  Table  A-5. 


Table  A-l.  ZONE  DEFINITION  ROUTINES 


Subroutine 

Use 

INIT 

Initialize  MMCLIB  variables 

SERINF 

Define  functions  and  material 
properties  to  be  used  In  zone 

MAPINF 

Define  mapping  function  In  zone 

CON INF 

Define  unknown  not  contained  In 
functions 

fir*  * 

15  j jfype 


Table  A-2.  NAPPING  FUNCTION  ROUTINES 


Subroutine 

Mapping  Function  Calculated 

Table  A-3. 

FUNCTION  COLLOCATION  ROUTINES 

NAPB 

-(c) 

HAP8AR 

Subroutine 

Function  Collocated 

MAPBR 

-tc) 

PHI 

of  (z) 

HAPBR1 

St  1/e) 

PHIB 

of  (D 

MAPB1 

-0/D 

PH1BP 

oTHz) 

MAPPBR 

w'(c) 

PHIBR 

at  (z) 

NAPPB1 

5T*(c) 

PHIB1 

of  (l/D 

MAPPB2 

rd/c) 

PHIB1P 

of'(l/z) 

NAPPP1 

-'■(e) 

PHIP 

of'(z) 

NAPP1 

-He) 

PHIPB 

of'(Z) 

MAPI 

-(c) 

PHIPP 

of(z) 

NAPB 12 

l/T,  u(l/e),  u2(-i(l/c)) 

PHIPPS 

of' ‘(z) 

NAPB21 

Z,  -(c).  -2(«l(c)) 

PHIBP1 

of'(D 

NAPPR 

1.  -He).  — i(— i(e) ) -He) 

PHIP1B 

of'(l/D 

MAPPRB 

1,  1,  —2  ( — 1 ( C ) ) 

CONSET 

unknown  not  In 

MAPPRP 

l.  l.  -H-He)) 

♦z  °r  *Z 

MAPP12 

MAPP21 

MAPPER 

I.  -'(!/?).  -H-lO/e))  «{(I/c) 

1«  -He).  -H-He))  -1(c) 

Z,  -(c).  -2(»l(c)) 

a Is  a complex  constant  and  f(z)  Is  a 
complex  function  (usually  either  4z(z), 
♦z(z),  4(c).  or  *(c)). 

z,  --He).  -iH-zHc)) 

Table  A-5.  RESULT  COMPUTATION  ROUTINES 


Table  A-4.  SPECIALIZED  EQUATION 
GENERATION  ROUTINES 


Subroutine 

Use 

FDCOL 

Collocation  of  displacements 
and  forces 

STCOL 

Collocation  of  stresses 

STITCH 

Equate  forces  and  displacements 
at  zone  boundaries 

COWIN 

Generation  of  equations  for  use 
In  combined  finite  element  - 
W1C  analysis 

Subroutine 

Use 

EVAL 

Evaluate  equation  with  solution 

FORES 

Calculate  forces  and  displacements 
with  solution 

STRRES 

Calculate  stresses  with  solutions 

PRTCOF 

Print  coefficients  solved  for 

PRTCON 

Print  unknown  found  which  Is  not 

In  functions 

PRNTMP 

Print  mapping  parameters 

PRNTMT 

Print  material  properties 

PRNTSR 

Print  series  Information 

PRTZON 

Print  all  Information  In  zone 

K1K2 

Calculation  of  stress  Intensities 
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APPENDIX  B.  SOURCE  LISTINGS  FOR  MMCLIB 


FUNCTION  BRANCH  (TH.THBR) 


FUNCTION  TO  FIND  ANGLE  TH*2.*N*PI  SUCH  THAT 
THBR  < ANGLE  < THBR+2.*PI 


COMMON  /PI FAC/  PI .HLFP1 , TWO PI 
TH1 *TH 

THP.THBR+TWOPI 

10  IF  (TH1 .GE.THBR)  GO  TO  20 
TH1.THU-TW0PI 
GO  TO  10 

20  IF  (TH1.LE.THP)  GO  TO  30 
TH1.TH1-TW0PI 
GO  TO  20 
30  BRANCH*TH1 
RETURN 
END 


SUBROUT ! NE  B4PHDP( XP . YP , I Z , CS2A , SN2A , I REF , A . B) 


SUBROUTINE  TO  CALCULATE  FUNCTION  MULTIPLYING 
SECOND  DERIVATIVE  OF  PHI  IN  STRESS  EQUATIONS 


CALL  MAPPR(XP,YP,XPR, YPR, IZ) 

XBR  = 0 . 

YBR=0. 

CALL  MAPPER  ( XP,YP,XZBR.,  YZBR.IZ.1  ) 
YZ3R=-YZBR 

GO  TO  (40, 10,20) , IREF 
10  CALL  MAPB21 (XP, YP.XBR, YBR, IZ) 

GO  TO  30 

20  CALL  MAPB12(XP, YP.XBR, YBR, IZ) 

30  YBR=-Y3r. 

40  XZBR=X2R-X2ER 
YZ8R«YBR-'ZBR 
XT*CS2A>XZBR-SN2A*YZBR 
YT=CS2A» YZBR+SN2A«XZBR 
XT1 =XPR*XPR-YPR*YPR 
YT1=2.*XPR*YPR 
DENOM= 1 ./(XT1*XT1+YT1*YT1 ) 

A = DEN0M*(XT*XT1+YT*YT1  ) 
B=DEN0M*(XT1*YT-YT1-XT) 

RETURN 

END 


SUBROUTINE  84PHIP  ( XP , YP , IZ , IOPT ,C2 , A , B) 


SUBROUTINE  TO  CALCULATE  FUNCTION  MULTIPLYING  COMPLEX 
CONJUGATE  OF  DERIVATIVE  OF  PHI  IN  FORCE , DISPLACEMENT 
EQUATIONS 
IOPT 

1 - NO  REFLECTION 

2 - X-AXIS  REFLECTION 

3 - UNIT  CIRCLE  REFLECTION 
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COMMON  /MAP/  IMAP( 5 , 2) , PARAM(  5 ,2, 5) 

CALL  MAPPER  (XP.YP.2X.2Y, IZ. 1 ) 

GO  TO  (40,10,20,60).  IOPT 
10  CALL  MAPB21  (XP.YP.ZXB.ZYB.1Z) 

GO  TO  30 

20  CALL  MAPB12  ( XP . YP . ZXB, ZYB. XZ) 

30  ZX«ZX-2XB 
ZY«ZY-ZYB 
40  XPRl » 1 . 

YPRl «0- 

IF  ( IMAP( IZ, 1 ) .EO.O)  GO  TO  50 
CALL  MAPPBR  (XP. YP, IZ. 1 . XPRl , YPRl ) 

50  CALL  MAPPRB  (XP,YP.XPR2,YPR2,IZ) 

DEN0M*C2/( ( XPR 1 *XPR1 +YPR1 * YPRl ) • ( XPR2*XPR2V*PR2*YPR2) ) 
2XB.XPR1 »XPR2-YPR1«YPR2 
ZYB»-XPR1*YPR2-XPR2*YPR1 
A*(ZX«ZXB-ZY-ZYB)*OENOM 
B=(ZY*2XB+ZX«ZYB) »DENOM 
RETURN 
60  A*C2*ZX 
B*C2*ZY 
RETURN 
END 


SUBROUTINE  84PHP(XP. YP.CS2A.SN2A. IZ, IREF.A.B) 


SUBROUTINE  TO  CALCULATE  FUNCTION  MULTIPLYING 
DERIVATIVE  OF  PHI  IN  STRESS  EQUATIONS 


I 

I 


COMMON  /MAP/  IMAP ( 5.2), PARAM( 5 ,2,5) 
A*0 . 

B=0 . 

if(imap(iz,i).ne.o)go  to  20 

IF(  IREF. EQ. 1 )G0  TO  80 
X2=1  . 

Y2  = 0. 

IF(IREF.EQ.2)G0  TO  10 

XT*XP*XP-YP*YP 

YT=2.*XP*YP 

DEN0M.1 ./(XT*XT+YT*YT) 

X2*-XT  *DENOM 
Y2=  YT*DENOM 
10  A=X2*CS2A-Y2*SN2A 
B=X2*SN2A+Y2*CS2A 
GO  TO  80 

20  CALL  MAPPER(XP,YP,X2,Y2,IZ,1) 

Y2—Y2 

CALL  MAPPR(XP, YP.OERX.DERY, IZ) 
XT=DERX»(DERX»DERX-3.»DERY*DERY) 
YT»OERY* (3. •DERX*DERX-DERY«OERY) 
DEN0M.1 ./(XT*XT+YT»YT) 

DX>xT  »DEN0M 
DY«-YT«DENOM 

CALL  MAPPP1(XP,YP.IZ,1 ,X1 ,Y1 ) 
IF(1MAP(1Z.2).EQ.O)GO  TO  30 
XT*X1 
VT«  Y 1 

CALL  MAPPRP( XP .YP.XT1 ,YT1  , I Z ) 

X1*XT»XT1-YT*YT1 

Y1*XT»YT1+YT*XTl 

CALL  MAPPl (XP.VP.IZ.1 ,XT1  ,YT1  ) 

XT«XT1»XTt-YT1*YT1 

YT«2 . »XT 1 • YT1 

CALL  MAPI (XP.YP.IZ.1, WX , WY ) 
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o o o o o 


CALL  MAPPP1 (WX.WY, IZ.2.XT1 ,Y11 ) 
XI *X1+XT*XT1-YT»YT1 
Y1 *Y1+XT*YT1+YT*XT1 
30  X3*0. 

Y3*0. 

GO  TO  (70,40,50) , IREF 
40  CALL  MAPB21 (XP.YP.XT.YT.IZ) 
YT*-YT 

CALL  MAPP21 (XP, YP.XT1 ,YT1 , IZ) 

GO  TO  60 

50  CALL  MAP PI 2(XP,YP,XT,YT,IZ) 
WX=XP*XP-YP«YP 
WY*2. *XP»YP 
DENOM* 1 ,/(WX*WX+WY*WY) 
XT1=-0EN0M*(WX*XT+WY*YT) 

YT 1 s-DENOM* ( WX* YT-WY » XT  ) 

CALL  MAPB12(XP,YP.XT,YT,IZ) 

YT  = -YT 
60  X2=X2-XT 
Y2=Y2-YT 

CALL  MAPPR (XP,YP,XT ,YT,IZ) 

X3  = XT  *XT 1 — YT* YT 1 
Y3=XT*YT1+YT»XT1 
70  XT=X2*XI-Y2*Y1 +X3 
YT=X2*Y1+Y2*X1+Y3 
XT1=DX«XT-DY«YT 
YT1=0X*YT+DY*XT 
A=CS2A*XT1-SN2A»YT1 
B=CS2A*YT1+SN2A*XT1 
80  RETURN 
END 


SUBROUTINE  COMBIN  ( NFI LE, IZ . IREF , IPRT , IERR) 

routine  to  write  equations  for 

COMBINED  FINITE  ELEMENT  - MMC 
ANALYSIS 


COMMON  /EON/  EQU( 3 , 1 00 ) , I BEG I (10),IENOI(10), EQUTOT (2,175), 

* SCALE ( 175) , NN , NUNKS , NUNK , NWRT , IUNI T 
COMMON  /MAT/  G(5),ETAM(5) 

COMMON  /MMCFE/  X ( 1 50 ) , Y ( 1 50 ) , FDRHS( 300 ) , EQNS( 1 20 , 1 76) , 

* KUV( 120,2) 

DIMENSION  RHS ( 1 ) 

REAL  KUV 

EQUIVALENCE  ( RHS( 1 ) , FQRHS( 1 ) ) 

OATA  NEQ  /1 00/ 

I CNT  = ICNT  + 1 

IF  (IERR.NE.O. AND. ICNT. EQ.1)  GO  TO  70 
DO  10  1*1 , NEQ 

DO  10  J*1 , 176 

10  EQNS( I , J ) *0 . 

REWIND  NFILE 

READ  (NFILE)  N 1 , N2 , I DOFST 

IF  (IPRT.EQ.2.AND.IERR.EQ.0)  PRINT  120.N1.N2 

NOOST  * ( I DOFST  + 1 )/2 

NPT.N2/2 

N00LST=NPT+N0DST-1 

READ  (NFILE)  ( X ( I ) , Y ( I ) , I * 1 , NPT ) 

READ  (NFILE)  ( RHS< I ) . I * 1 , N I ) 

DO  50  N*1 .NPT 

CALL  MAPPER(XP, YP,X(N) ,Y(N) , IZ,2) 

IF  (IPRT. NE. 2. OR. IERR.NE.O)  GO  TO  20 
N0DE=N+N00ST-1 

PRINT  130,  NODE ,X(N) ,Y(N) ,XP,YP 
20  CALL  ZERO  ( 0 . 0 , 1 . 2 . 1 , 0 ) 

CALL  PTCOL  (XP,YP,IZ,2,IREF,1.,1,2,1,2) 


oooooooooo 


IF  (IREF.EQ.4)  call  MLTOSP  (XP.VP.IZ) 

DO  30  J*1 .2 

30  REAO  (NFUE)  (KUV(  I,J),I*1,N1) 

00  40  1*1 ,N1 

DO  40  J«1 , NUNK 

40  E0NS(I.J)*E0NS(I,J)+KUV(I,1)*E0UT0T(1 , J)+KUV ( l ,2) *EQUT0T (2 ,U) 
50  CONTINUE 

TW0G*2 • *G( IZ) 

DO  60  1*1, N1 

60  EQNS(I.NN)*RHS(I)*TWOG 

IF  (IPRT.GT.O.AND.IERR.EO.O)  PRINT  140,  N1 .NOOST .NODLST 
70  N1HlF.N1/2 

DO  100  N.1.N1HLF 

I E2*2*N 

IE1.IE2-1 

DO  80  I = 1 , NN 

EOUTOT(1,I)*EONS(IE1,I) 

80  EOUTOTI 2 , I ) >EQNS( I E 2 . I ) 

IF  ( IERR.NE.0)  GO  TO  90 
CALL  SAVE  (1) 

GO  TO  100 

90  CALL  EVAL  (1 ,2,2.XER,YER) 

PRINT  110,  XER.YER 

100  CONTINUE 
RETURN 

110  FORMAT  ( 1 H .2F10.4) 

120  FORMAT  1////5H  N1  * , 1 5 , 1 0X , 4HN2  *,I5 

* /////36H  EQUATIONS  FOR  COMBINED  METHOD  TAKEN, 

1 28H  FROM  FOLLOWING  NODAL  POINT S//2X , 4HNOOE , 1 1 X , 

2 1HX.11X.1HY) 

130  FORMAT  (1H  .I5.4E12.4) 

140  FORMAT  (////I6.33H  EQUATIONS  GENERATED  ON  INTERFACE/ 

* 9H  AT  N00ES.I5.5H  THRU, 15) 


SUBROUTINE  CONINF( ICONST, IRLIM) 


f 


SUBROUTINE  TO  SET  UP  FOR  EXTRA  UNKNOWNS  SUCH  AS  FORCE  CONSTAN 
ICONST  - NUMBER  OF  EXTRA  UNKNOWN 
IRLIM  - TWO  DIGIT  NUMBER 

01  — UNKNOWN  HAS  ONLY  IMAG  PART 

10  — UNKNOWN  HAS  ONLY  REAL  PART 

11  — UNKNOWN  HAS  BOTH  REAL,  IMAG  PARTS 


DIMENSION  I IMR L ( 2 ) 

EQUIVALENCE  ( I RL , I IMRL ( 2 ) ) , ( I IM , I IMR L ( 1 ) ) 

COMMON  /CONST/  ICON (5,2),ISTRT(5) 

COMMON  /EQN/  EQU ( 3 , 1 00 ) , I BEG I ( 1 0) , I END I (10), EQUTOT (2,175), 
• SCA  LE ( 1 75 ) , NN , NUNKS , NUNK , NWRT , I UN IT 

CALL  DECODE( IRLIM, 2, IJMRL.NDUM) 

IF(IRL.GT.0)IRL«1 
IF(IIM.GT.0)IIM»1 
NADD* IR  L+ 1 IM 
IF(NADD.LE.O)GO  TO  10 
ICON ( ICONST ,1 ) • I RL 
ICON ( ICONST , 2 ) * I IM 
ISTRT( ICONST ) *NUNK+1 
NUNK *NUNK+N ADD 
NNsNUNK+1 


10  RETURN 
END 


ci n n oooono  nnnnnnnnnn 
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SUBROUTINE  CON SET ( ICNST , XR 1 , XR2 , XI 1 , XI2 ) 


SUBROUTINE  TO  SET  VALUES  FOR  COLLOCATION  OF  AOOEO  UNKNOWNS 
ICNST  - ADDED  UNKNOWN  NUMBER 

X91  - VALUE  FOR  real  EQUATION,  REAL  PART  OF  UNKNOWN 
XR2  - VALUE  FOR  REAL  EOUATION,  IMAG  PART  OF  UNKNOWN 
XII  - VALUE  FOR  IMAG  EQUATION,  REAL  PART  OF  UNKNOWN 
XI2  - VALUE  FOR  IMAG  EQUATION,  IMAG  PART  OF  UNKNOWN 


COMMON  /CONST/  I CON (5, 2) , ISTRT(5) 

COMMON  /EQN/  EQU( 3, 1 00 ) , I BEG I (10), I END I ( 10), EQUTOT(2, 175) , 
* SCALE( 175) .NN.NUNKS.NUNK.NWRT, IUNIT 

11*1 STRT ( ICNST ) 

!F( ICON( ICNST, 1 ) .LE.OJGO  TO  10 
EQUTQT ( 1 . II )*XR1 
EQUTOT ( 2,I1)*XI1 
11*11+1 

10  IF(ICON(ICNST,2).LE.O)GO  TO  20 
EQUTOT ( 1 , 1 1 )*XR2 
EQUTOT (2,I1)*XI2 
20  RETURN 
END 


subroutine  crsmlt 


SUBROUTINE  TO  READ  EQUATIONS  FROM  TAPE  AN0  PERFORM 
SCALING,  CROSS  MULTIPLICATION 


COMMON  /EQN/  EQU (3,100),IBEGI(10),IENDI(10), EQUTOT (2,175), 
* SCALE( 175) .NN.NUNKS.NUNK.NWRT, IUNIT 

COMMON  /SOLN/  A( 1 540 1 ) . 1 01 AG( 1 76 ) , PHI ( 1 76) 

REWIND  IUNIT 
DO  40  1*1 , NWRT 

READ  (IUNIT)  ( ( EQUTOT  («J,K)  ,K*1  ,NN)  ,d*1 ,2) 

DO  10  d*1 ,2 

DO  10  K= 1 , NUNK 

10  EQUTOT(U,K )= EQU TOT (J,K)*SCAL£(K) 

DO  30  L* 1 ,2 
IND*0 

DO  20  J=1 , NUNK 

00  20  K* J , NUNK 
I ND* I ND+ 1 

20  A ( IND)*A(INO)+EQUTOT(L,J)*EQUTOT(L.K) 

DO  30  U* 1 , NUNK 

30  PHI (J)*PHI(J ) + EQUTOT ( L , J) * EQUTOT ( L , NN ) 

40  CONTINUE 
RETURN 
END 


SUBROUTINE  DECOOE( lOPT .MAX . IDIGIT ,NUM) 
ROUTINE  TO  SEPARATE  DIGITS  FFROM 
INTEGER  CONTAINING  STRING  OF  DIGITS 


DIMENSION  IDIGIT ( 1 ) 

NUM*0 

ITSTRsI 

DO  10  1*1 .MAX 
10  1 0 IGI T ( I ) *0 
DO  20  1*1 .MAX 
IF(ITSTR.GT.IOPT)GO  TO  30 
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NUMxNUM+1 
IDEN* ITSTR 
I TSTR*1 0*« I 

20  I DIGIT ( I )«MOO( IOPT , ITSTR) /IDEM 
30  RETURN 
END 


SUBROUTINE  DOF(NDOF) 

ROUTINE  TO  INTERROGATE  MMCLIB  DATA  BLOCKS 
FOR  NUMBER  OF  UNKNOWNS 

COMMON  /EON/  EQU( 3 ,100), IBEGI ( 1 0) , I ENDI ( 10), E0UT0T(2 . 175) 
SCALE ( 1 75 ) , NN .NUNKS , NUNK.NWRT, I UN IT 

ND0F*NUNK 

RETURN 

END 


SUBROUTINE  EVAL  ( IE  1 , 1 E2 , IOPT , VAL1 , VAL2 ) 


SUBROUTINE  TO  EVALUATE  EQUATIONS  ONCE  SOLUTION  IS  KNOWN 


COMMON  /EON/  EQU( 3 , 1 00 ) , I BEG I (10),IENDI(10), EOUTOT (2,175) 
• SCALE ( 175) , NN, NUNKS, NUNK.NWRT, IUNIT 

COMMON  /SOLN/  A ( 1 540 1 ) , I DI AG( 1 76 ) , PHI ( 1 76) 

DIMENSION  VALUE ( 2 ) 

DO  10  IE=IE1 . IE2 
VALUE(IE)*0. 

IF  (I0PT.E0.2)  VALUE(IE)*-EQUTOT(IE,NN) 

DO  10  J = 1 , NUNK 

10  VALUE(IE)=VALUE(IE)+PHl(U)*EQUTOT(IE,d) 

VAL 1 * VA  LUE ( t ) 

VAL2*VALUE(2) 

RETURN 

END 


SUBROUTINE  FDRES  ( IZ , NPT , IOPT , IREF ,X , Y ,XP,TP) 


SUBROUTINE  TO  PRINT  FORCES,  DISPLACEMENTS  FROM  SOLUTION 


IOPT 


1 FOR  FORCES 

2 FOR  DISPLACEMENTS 

1 FOR  NO  REFLECTIONS  IN  ZONE 

2 FOR  X-AXIS  REFLECTION 

3 FOR  UNIT  CIRCLE  REFLECTION 


IREF 


DIMENSION  X ( 1 ) . Y( 1 ) , XP(1),  YP( 1 ) 
COMMON  /MAT/  G(5).ETAM(5) 

A*  1 . 

B*0 . 

GO  TO  (10,20),  IOPT 
10  PRINT  50 
GO  TO  30 
20  PRINT  60 
30  DC  40  1*1 , NPT 

CALL  ZERO  (0,0,1 ,2, 1.0) 

CALL  PTCOL  (XP(I),YP(I),IZ. IOPT, IREF, 1 


r 


IF  (IREF.EQ.4)  CALL  MAPPR  { XP ( I ) , YP ( I ) , A , B , IZ) 

FACT  * 1 . / ( A*A+B»B) 

IF  (IOPT.E0.2)  FACT=FACT/(2.*G(IZ)) 

CALL  EVAL  ( 1 ,2, 1 ,V1 , V2) 

VALX=(V1  » A-V2 • B ) » FACT 
VALY*(V1 »B+V2»A) »FACT 

40  PRINT  70,  VALX1VALYtX(I),Y(I)tXP(I),YP(I) 

RETURN 

50  FORMAT  (21HOEVALUATION  OF  FORCES  / 9X , 2HF 1 . 1 1 X . 2HF2 , 1 2X , 1HX , 1 2X , 1 H 
1 Y 8X , 8HX (PARAM) ,5X,8HY(PARAM)  ) 

60  FORMAT  ( 28H0EVALUAT ION  OF  DISPLACEMENTS  / 1 OX , 1 HU , 1 2X , 1HV , 1 2X , 1HX . 

1 1 2X , 1 H Y , 8X , 8HX ( PARAM) , 5X , 8HY ( PARAM)  ) 

70  FORMAT  (1X.6E13.5) 

END 


SUBROUTINE  FDCOL  (XP.YP.N, 12, IOPT, IREF, RHSt , RHS2, I ERR) 

C 

C 

C SUBROUTINE  TO  COMPUTE  FORCE,  DISPLACEMENT  EQUATIONS 

C MEANING  OF  IREF 

C 1 - NO  REFLECTIONS 

C 2 - X-AXIS  REFLECTION 

C 3 - UNIT  CIRCLE  REFLECTION 

C 
C 

DIMENSION  XP(1),  YP ( 1 ) , RHSl(t),  RHS2(1) 

COMMON  /EON/  EQU( 3 , 1 00 ) , I BEG I (t0),IENDI(10), E0UT0T(2 , 175), 
• ' SCALE ( 175) ,NN, NUNKS , NUNK , NWRT , I UNIT 

PRINT  40 
A*  1 . 

B*0. 

DO  20  I«1  , N 

CALL  ZERO  (0,0,1 ,2,1 .0) 

CALL  PTCOL  ( XP ( I ) , YP(I), IZ, IOPT, IREF, 1., 1,2, 1,2) 

I F ( IREF.EQ.4)CALL  MAPPR(XP( I ) , YP ( I ) , A , B, IZ ) 

EOUTOT ( 1 , NN ) a A*RHS 1 ( I )+B*RHS2( I ) 

EQUTOT ( 2 , NN ) * A*  RHS2 ( I )-B*RHS1 ( I ) 

IF  (IERR.E0.1)  GO  TO  10 
CALL  SAVE  (1) 

GO  TO  20 

10  CALL  EVAL  (1 ,2,2,XER,YER) 

PRINT  30,  XER,YER,XP(I),YP(I) 

20  CONTINUE 
RETURN 

30  FORMAT  ( 1 H , 2F 1 0 . 4 , 2E1 0 . 4 ) 

40  FORMAT( 1 H ) 

END 


SUBROUTINE  FNOMAT 

c 

C 

C SUBROUTINE  TO  SCALE,  CROSS-MULTI  PLY , ANO  SOLVE 

C FOR  COEFFICIENTS 

C 

c 

COMMON  /eon/  EQU( 3 , 1 00 ) , I BEG  I (10),IENDI(10), EQUTOT (2 ,175), 
« SCALE ( 175) , NN, NUNKS, NUNK, NWRT, IUNIT 

COMMON  /SOLN/  A ( 1 540 1 ) , I DI AG( 1 76 ) , PH  I ( 1 76 ) 

CALL  SLVINT 
CALL  SCLFIX 
CALL  CRSKLT 
CALL  SOLVE 
DO  10  1=1 , NUNK 

10  PHI(I)*PHI(I)*SCALE(I) 

RETURN 

END 
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SUBROUTINE  INIT 


SUBROUTINE  TO  INITIALIZE  ALL  COMMON  VARIABLES 


COMMON  /CONST/ 
COMMON  /EON/ 

COMMON  /MAP/ 
COMMON  /MAT/ 
COMMON  /PI  FAC/ 
COMMON  /POWERS/ 
COMMON  /SERIES/ 


COMMON  /SOLN/ 
PI. 3. 1^159265 
HLFPI*1 .57079633 
TW0PI.6. 28318531 
IUNIT.7 


I CON (5, 2) , ISTRT(S) 

EQU( 3 , 100 ) , I BEG I ( 10) , I END I ( 10), E0UT0T(2  >175) 
SC A L E ( 1 75  > , NN . NUNKS , NUNK , NWR  T , I UN I T 
IMAP ( 5 , 2) , PaRAM( 5 ,2,5) 

0(5) , E T AM ( 5 ) 

PI.HLFPI.TWOPI 

R( 1 00) , TH( 1 00) , ISHFT 

ISTYP( 10) . IZONES , XNOT (10) , VNOT( 10) , NNEG( 10) , 
NPOS (10),IZRO(10), IOORL( 10). 
IODIM(10),IEVRL(10).IEVIM(10) 

A( 1 5401 ),I0IAG( 178), PHI (176) 


ishft.o 
00  10  1*1,10 
IBEGI ( I ) *0 
1ENDI ( I ) =0 
NNEG( I ) *0 
NPOS(I)*0 
IZRO( I ) *0 
IODRL( I ) -0 
IOOIM( I ) *0 
IEVRL(I)«0 
IEVIM(I).0 
XNOT ( I ) =0. 

10  YNOT ( I ) *0, 

DO  20  1*1,5 

1STRT ( I ) *0 
G( 1 ) *0 . 

ET  AM( I ) *0 . 

00  20  J*1  ,2 

IMAP( I . J)*0 
ICON(I.J)*0 
00  20  K*1 ,5 

20  PARAM(I.J,K)«0. 
REWIND  I UN IT 


IZONES*0 

NWRT*0 


NUNK*0 

nunks*o 

NN*0 

00  30  1*1 ,175 

30  SCALE{I)*0. 
RETURN 
END 


SUBROUTINE  K1K2(XP, YP, A,B,K1 ,K2,tZ) 

ROUTINE  TO  COMPUTE  STRESS 
INTENSITY  FACTORS 

REAL  K1 ,K2 

CALL  ZERO( 1 ,2 , 1 , 2 ,3 , 0) 

CALL  PHIPB(XP,YP,A.-B, 1Z) 

CALL  PLACEUZ.1  ,2.1  .2, 1 .) 

CALL  EVAL(1,2,1,K1,K2) 

RETURN 

END 
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SUBROUTINE  MAP  (XP, YP, PI , P2, P3, P4, P5.M.X, Y) 


SUBROUTINE  TO  CALCULATE  MAPPING  FUNCTIONS 


10  CSPHI«C0S(P3) 

SNPH I *S 1 N( P3 ) 

RS0*XP*  XP+YP* YP 
RS0S-RS0+1 . 

RS0D-RS0-1 . 

SUMl * (PI *RSOS+P2*RSOO) *XP 
SUM2* ( P 1 •RS00+P2*RS0S) *YP 
DEN0M.1 ./(2.*RSQ) 

X^P4+(CSPHI*SUM1-SNPHI*SUM2)*0EN0M 
Y • P5+ ( SN PH I • SUM 1 ♦ CS PH I * SUM2 ) • OENOM 
GO  TO  30 


MAPPING  *2 


20  CSPHl*C0S( P3) 

SNPhI*SIN(P3) 

snhy*sinh(yp) 

CSHY  *C05H( YP ) 

SUMl*SIN(XP)*( P2*SNHY+P1 *CSHY) 
SUM2-C0S ( XP ) ♦ ( P2*CSHY+P1 « SNHT ) 
X*P4+CSPHI*SUM1-SNPHI*SUM2 
Y«P5+SNPHI •SUM1+CSPHI*SUM2 
30  RETURN 
END 


SUBROUTINE  MAPB  (XP.YP.IZ. I12.X, Y) 

CALL  PARA  (IZ,I12.M,P1 ,P2,P3,P4,P5) 
CALL  MAP  (XP.-YP.P1 ,P2,P3.P4,P5,M,X,V) 
RETURN 
END 


SUBROUTINE  MAPBAR  (XP.YP.IZ, 112, X,Y) 

CALL  PARA  (IZ.I12.M,P1,P2.P3.P4.P5) 

CALL  MAP  (XP,YP,P1.P2,P3.P4,P5.M,X,Y) 

Y— Y 

RETURN 

END 


SUBROUTINE  MAPBR  (XP.YP.IZ, 112, X.Y) 
CALL  PARA  (IZ,I12,M,P1 ,P2,P3,P4,PS) 

GO  TO  (10,20),  M 
10  BB»P2 
GO  TO  30 
20  BB*-P2 

30  CALL  MAP  (XP.YP.P1 ,BB,P3,P4,-P5,M,X,Y) 
RETURN 
END 


SUBROUTINE  MAPBR1  ( XP , YP , I Z , 1 1 2 , X , Y ) 
CALL  MAPB t (XP.YP.IZ, 112, X.Y) 

Y*“  Y 

RETURN 

END 


uuuuu 


SUBROUTINE  MAPB1  (X?,YP,IZ.I12.X,Y) 
CALL  PARA  (IZ.112.M.P1.P2.P3.P4.P5) 

R1 * 1 ,/(XP*XP*YP.YP) 

XT»XP»R1 

YT*YP«R1 

CALL  MAP  (XT.YT.P1 ,P2,P3,P4,PS.M,X.Y) 

RETURN 

ENO 


SUBROUTINE  MAPB12  ( XP . YP, X . V , IZ) 

COMMON  /MAP/  IMAP(S,2) . PARAM( 5, 2, 5) 
R1 =1 ,/(XP*XP+YP»YP) 

X»XP«R1 
Y* YP»R1 

IF  (IMAP(IZ.I).EO.O)  RETURN 
CALL  MAPB1  (XP.YP.IZ.1 .X.Y) 

IF  (IKAP(IZ,2).EQ.0)  RETURN 

XX«X 

YY*Y 

CALL  MAPI  (XX , YY , IZ . 2 , X, Y) 

RETURN 

ENO 


SUBROUTINE  MAPB21  ( XP , YP , X , Y , IZ ) 

COMMON  /MAP/  IMAP(5.2) . PARAM( 5 , 2 . 5) 

X*XP 

Y=-YP 


IF  (IMAP(IZ.I).EO.O)  RETURN 
CALL  KAPB  (XP.YP.IZ.1, X.Y) 

IF  ( IMAP ( IZ,2) . EO. 0 ) RETURN 

XX«X 

YY*Y 

CALL  MAPI  (XX , YY , IZ . 2 , X, Y) 

RETURN 

ENO 


SUBROUTINE  MAPINF  (IZ.M1.P11 ,P12,P13,P14,P15,M2,P21,P22 
1 P23.P24.P25) 


SUBROUTINE  TO  SET  MAPPING  FUNCTION  VARIABLES 


COMMON  /MAP/  IMAP(5 ,2) , PAR AM (5, 2, 5) 
DATA  ANGCON  / . 1 7453292E-1 / 

IF  (M1.LE.O)  RETURN 
PARAM( IZ. 1 , 1 )«P1 1 
PARAM( IZ. 1 ,2)*P12 
P ARAM( I Z , 1 , 3) *P1 3* ANGCON 
PARAM( IZ, 1 , 4 ) » PI  4 
PAR AM( I Z . 1 ,5)»P15 
IMAP( IZ, 1 )«M1 
IF  (M2.LE.O)  RETURN 

par am ( iz,2.i)*P2i 

PARAM( IZ.2,2)>P22 

PARAM( I Z , 2 , 3 ) *P23* ANGCON 

PARAM( IZ , 2 , 4 ) * P24 

PAR AM( IZ.2,5)«P25 

IMAP( IZ , 2) «M2 

RETURN 

ENO 
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SUBROUTINE  MAPP  (XP.YP.P1 ,P2,P3.P4,P5,M,X0.VD) 


routine  to  compute  first  derivative 
of  mapping  function 

COMMON  /MAP/  IMAP(5,2) ,PARAM(5,2,5) 

COMMON  /P I FAC/  Pl.HLFPI.TWOPI 
GO  TO  (10,20),  M 

MAPPING  #1 

10  BMA.P2-P1 

R4»XP-XP+YP»YP 
R4aR4«R4 
R41.1 ./R4 

FACT1«.5*(BMA»(XP*XP-VP*YP)+R4*(P1+P2) ) 

FACT2*BMA.XP»YP 

CS«C0S(P3) 

SN>SIN( P3 ) 

XD*R41 • ( FACT1  •CS+FACT2*SN) 
YD*R41*(FACT1»SN-FACT2*CS) 

RETURN 

20  CALL  MAP  (XP.YP.P2.P1 .P3-HLFPI.0. ,0. .M.XO.YO) 
RETURN 

ENO 


SUBROUTINE  MAPINV  (X.Y.P1 ,P2.P3,P4,P5,M.XP,VP) 


SUBROUTINE  TO  CALCULATE  MAPPING  FUNCTION  INVERSES 


COMMON  /PIFAC/  PI.HLFPI.TWOPI 
GO  TO  (10,20).  M 

MAPPING  #1 

10  CALL  QUADF  (PI , P2 ,P3 , P4 ,P5 , X , Y , 1 ,XP , YP) 

RETURN 

MAPPING  #2 

20  IMAPR«2 

30  CALL  QUADF  ( P2 ,-Pi , P3+HLFPI , P4, P5 , X , Y , IMAPR , XX, YY) 
XP.XX 
VP«YY 

IF  (P1.EQ.P2)  GO  TO  50 
XP*ATAN2( YY.XX) 

YP«-ALOG(SQRT(XX*XX-FYV*YY)  ) 

IF  (VP. GT.-. 0001 ) GO  TO  SO 
IF  ( IMAPR. NE. -2)  GO  TO  40 
WRITE  (6,60)  M.X.Y 
STOP 

40  IMAPR<-2 
GO  TO  30 

50  RETURN 

60  FORMAT  (17H  ERROR  IN  INVERSE  .I3.2E12.S) 

ENO 


SUBROUTINE  MAPPB1  ( XP, YP, IZ , 1 12 , XD , YD) 

CALL  PARA  (IZ, I12.M.P1 ,P2,P3,P4,P5) 

GO  TO  (10,20),  M 

10  CALL  MAPP  ( XP,  YP.P1 ,P2,-P3,P4,P5.M,XD,YD) 
RETURN 

20  CALL  MAPP  (-XP ,-YP, PI , P2 , -P3 , P4 , PS ,M, XD, YD) 
RETURN 
END 
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SUBROUTINE  MAPPBR  (XP,YP,IZ,I12.XD,Y0) 
CALL  PARA  (IZ.H2.M.P1  .P2.P3.P4.P5) 

CALL  MAPP  (XP.YP.P1 .P2,P3.P4,P5.M,X0.YD) 

YD— YD 

RETURN 

END 


SUBROUTINE  MAPP82  ( XP , YP, IZ , I 12 .XO.YO) 
RR« 1 ./SQRT(XP«XP+YP»YP) 

XX>XP*RR 

YY»-YP*RR 

CALL  MAPPB1  (XX.YY.IZ.I12.XD.YD) 

RETURN 

END 


SUBROUTINE  MAPPOP(XP,YP,XD.YO.IZ) 

COMMON  /MAP/  IMAP(5 , 2) , PAR AM (5,2 ,5) 

XD*1  . 

YD«0. 

I F ( IMAP ( IZ.1 ).EQ.O.AND. IMAP( IZ. 2) . EQ.O) RETURN 

CALL  MAPPP1 (XP.YP.IZ.1 .ZMX2.ZMY2) 

IF(IMAP(IZ.2).NE.O)GO  TO  10 

XD«ZMX2 

Y0*ZMY2 

RETURN 

10  CALL  MAPP1  (XP.YP.IZ.I, ZMX 1 , ZMY1 ) 

ZMX 1 SQ= ZMX 1 » ZMX 1 -ZMY 1 * ZMY 1 

ZMY1S0=2.*ZMX1»ZMY1 

CALL  MAPB(XP,-YP, IZ.WXP.WY) 

CALL  VAPP1  (WX.WY.IZ.2.WMX1 , WMY 1 ) 

CALL  MAPPP1 ( WX , WY , I Z , 2 ■ WMX2 , WMY2 ) 

XD»  WMX2 • ZMX 1 SQ- WMY2  * ZMY  t SQ+WMK1 • ZMX2-WMY 1 *ZMY2 
YD>  WMX2 * ZMY  1 SQ+WMY2* ZMX  1 SQ+WMX 1 • ZMY2+WMY 1 *ZMX2 
RETURN 
END 


SUBROUTINE  MAPPP1  ( XP , YP. IZ . 1 1 2 . XD2 , YD2 ) 
CALL  PARA  (IZ.I12.M.P1 .P2.P3.P4.P5) 

CALL  MAPPP  ( XP , YP , P 1 .P2.P3.P4, P5 ,M, X02 , YD2) 

RETURN 

END 


SUBROUTINE  MAPPR(XP. YP.XD. YD, IZ) 

COMMON  /MAP/  IMAP(5,2) , PARAM( 5 , 2 , 5) 
XO*  1 . 

YD*0 . 

IF( IMA?( IZ, 1 ) .EO.OJGO  TO  10 
CALL  MAPP1 (XP.YP, IZ.1 ,XP1 ,YP1 ) 

CALL  KAPPRP(XP, YP.XP2.YP2, IZ) 
XD=XP1*XP2-YP1*YP2 
Y0*XP1»YP2+XP2*YP1 
10  RETURN 
END 


J 
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SUBROUTINE  MARRED  (XR.VR.X.V. 1Z, INV) 


SUBROUTINE  TO  CALCULATE  MARRING  FUNCTIONS  IN  ZONES 


COMMON  /MAR/  IMAR(5,2) , RARAM(S,2 .5) 

XX. XP 
YY«YP 
XXM.X 
YYM.Y 
IMP. 2 

IF  (IMAP(IZ.I).GT.O.ANO.  1MAR( IZ.2) .GT.O)  GO  TO  20 
IMP.  1 

IF  (IKAP(IZ,1 J.GT.O)  GO  TO  30 
GO  TO  (10,20).  INV 
10  X.XP 
Y.YP 
RETURN 
20  XP.X 
YP*  Y 
RETURN 

30  DO  60  I>1 .IMP 

II. I \ 

IF  ( INV. EQ. 2. AND. IMP.E0.2)  11*3-1 

CALL  PARA  (IZ.1I,M.P1,P2,R3,R4.RS) 

GO  TO  (40.50),  INV 

40  CALL  MAP  (XX.YY.P1 .P2,P3,P4,P5.M,X,V) 

XX* X 
YY»Y 

GO  TO  60 

50  CALL  MAP INV  (XXM.YYM.P1 ,R2,R3.P4,P5,M,XR.YR) 

XXM.XP 
YYM.YP 
60  CONTINUE 
RETURN 
END 


SUBROUTINE  MAPP1  ( XP , YP , IZ , 1 1 2. XD. YD) 
CALL  PARA  ( IZ , 1 1 2 ,M, Pi , P2 , P3« P4, R5) 

CALL  MAPP  (XP. YP.P1 ,P2,P3.P4,P5,M,XD,YD) 

RETURN 

END 


SUBROUTINE  MAPPP(XP . YR, PI , P2 , P3 , P4, P5.M, XD2. Y03) 

ROUTINE  TO  COMPUTE  SECOND  DERIVATIVE 
OF  MAPPING  FUNCTION 

GO  TO  (10,20).  M 

MAPPING  41 

10  R6.XP*XP+YP.YP 

R6.(P1-P2)/(R6«R6*R6) 

FACT1 »XP.XP*XP-3 .»XP*YP.YP 
FACT2.YP*  YP* YP~3 .*YP*XP*XP 
CS=C05( P3) 

SN»SIN( P3) 

XD2-R6* ( FACT1 .CS-FACT2.SN ) 

YD2-R6* ( FACT  1 *SN+FACT2»CS) 

RETURN 

MAPPING  42 
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ao  CALL  MAP  (XP,VP,P1.P2,P3.O.,0..M.XO2.VD2) 
XD2--XD2 
YD2--VD2 
RETURN 
END 


SUBROUTINE  MAPPRS  (XP, YP. XD, VO,  IZ) 
CALL  NAPPRP  (XP, VP.XD.VD, IZ) 

VO--YO 

RETURN 

END 


SUBROUTINE  MAPPRP  { XP, YP, XD, YD, IZ) 

COMMON  /MAP/  IMAP (5, 2) , PAR AM ( 5, 2, 5) 

XO-1  . 

Y0«0. 

IF  ( IMAP( IZ, 1 ) . EO.O . OR . IMAP( IZ,  2 ) . EQ . 0 ) RETURN 
CALL  MAPI  (XP.YP.IZ.I.XX.YY) 

CALL  MAPP1  (XX, YY , IZ ,2 ,XD, YD) 

RETURN 

END 


SUBROUTINE  MAPP1 2(XP. YP.XD. YO, IZ) 
DENOMxl . /( XP*XP+YP*YP) 

X1-XP-OENOM 

YIxyp.OENOM 

CALL  MAPPR(X1,Y1,XP.VD.1Z) 

YO—YD 

RETURN 

END 


SUBROUTINE  MAPP2I (XP,YP,XD,YD,IZ) 
CALL  MAPPR(XP.-YP,XD.YO,IZ) 

YD«-VO 

RETURN 

END 


SUBROUTINE  MAPI  (XP, YP, IZ, I 12, X, Y) 
CALL  >ARA  (IZ.I12,M.P1,P2.P3,P4,P5) 
CALL  MAP  (XP, YP, PI , P2, P3, PA,P5,M,X. V ) 
RETURN 
ENO 


SUBROUTINE  MLTDSP(XP.YP.IZ) 
CALL  MAPPR(XP.YP,A,B,IZ) 
AA*A/(A*A+B*B) 

BB*B/(A*A+B»B) 

CALL  R0TZ0N(3.IZ,AA,-BB,8B,AA) 

RETURN 

ENO 
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subroutine  nwrtfl(n) 


SUBROUTINE  TO  ALTER  VARIABLE  CONTAINING  NUMBER 
OF  EQUATIONS  WRITTEN  ON  CURRENT  SCRATCH  FILE 


COMMON  /EQN/  EQU( 3 , 1 00 ) , I BEG I ( 1 0 ) , I ENOI (10), EQUTOT (2.175) 
SCA  LE ( 1 75 ) , NN , NUNKS . NUNK , NWR  T , I UN I T 

NWRT  *N 
RETURN 
END 


FUNCTION  NWRTSV(IOUM) 

FUNCTION  TO  SAVE  NUMBER  OF  EOUATIONS  WRITTEN 
ON  CURRENT  SCRATCH  UNIT 

COMMON  /EQN/  EQU(3.100).IBEGI(10),IENDI(10),EQUTOT(2.17S) 
SCALE ( 1 75 ) , NN , NUNKS , NUNK , NWRT , I UNI T 

NWRTSV«NWRT 

RETURN 

END 


SUBROUTINE  para  (IZ. I12,M, PI ,P2, P3.PA.P5) 


SUBROUTINE  TO  FINO  MAPPING  FUNCTION  PARAMETERS  IN  ZONE 


COMMON  /MAP/  IMAP (5, 2) . PARAM( 5 , 2 , 5 ) 
M * IMAP  ( IZ. 1 12) 

P1*PARAM(IZ, 112,1) 

P2*PARAM( IZ, 1 12,2) 

P3*PARAM( IZ , 1 1 2 , 3) 

P4«PARAM( IZ, 112,4) 

P5*PARAM(IZ,I12,5) 

RETURN 

END 


SUBROUTINE  PHI  (XP.YP.A.B.IZ) 

routine  to  collocate  a function  in  a zone 


CALL  ZETN  (XP.YP.IZ) 
call  SERFIL  (lz,1..t 

RETURN 

END 


SUBROUTINE  PHIB  (XP.YP.A.B.IZ) 
CALL  PHI  (XP,-YP.A,B,IZ) 

RETURN 

END 


SUBROUTINE  PHI  BP  (XP.YP.A.B.IZ) 
CALL  ZETN1  (XP.YP.IZ) 

CALL  SERFIL  ( IZ  . 1 . ,-1 . . A , B) 

RETURN 

END 


oooo  noon 


SUBROUTINE  PHIBP1(XP,YP, A.B.IZ) 
CALL  PHIP(XP.-YP,A,B,1Z) 

RETURN 

END 


SUBROUTINE  PHI BR(XP , YP , A, B, 12) 
CALL  ZETN(XP, YP, IZ) 

CALL  SERFIL(IZ,-1.,-1..A,B) 

RETURN 

END 


SUBROUTINE  PHIB1  (XP.YP.A.B.IZ) 
R1«1 ./(XP.XP+YP.YP) 

CALL  PHI  (XP.R1 ,YP.R1 .A.B.IZ) 

RETURN 

END 


SUBROUTINE  PH IBlP(XP,YP, A.B.IZ) 
RR»  1 . /( XP. XP+YP. YP ) 

CALL  PHI  BP  (XP.RR.-YP.RR, A.B.IZ) 

RETURN 

END 


SUBROUTINE  PHI P ( XP . YP , A , B , 12) 

ROUTINE  TO  COLLOCATE  DERIVATIVE 
OF  A FUNCTION 

CALL  ZETN1 (XP.YP.IZ) 

CALL  SERFIL(I2,1.,1.,A,B) 

RETURN 

END 


SUBROUTINE  PHIPB(XP,YP, A.B.IZ) 
CALL  ZETN1 (XP.YP.IZ) 

CALL  SERFIL( IZ ,-1 . ,-1 . , A, B) 

RETURN 

END 


SUBROUTINE  PHIPP  (XP.YP.A.B.IZ) 

ROUTINE  TO  COLLOCATE  SECOND 
DERIVATIVE  OF  A FUNCTION  IN  a ZONE 

CALL  ZETN2  (XP.YP.IZ) 

CALL  SERFIL  ( I Z , 1 . , 1 . , A. B) 

RETURN 

ENO 
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SUBROUTINE  PHIPPB(XP,VP,A.B, 12) 
CALL  ZETN2(XP,YP,IZ) 

CALL  SERFIL(IZ,-1 .,-1. ,A.B) 

RETURN 

END 


subroutine  phipib(xp,yp.a.b, iz) 

RR»1 ,/(XP«XP+YP*YP) 

CALL  PHIP(XP*RR,YP*RR,A,B,IZ) 

RETURN 

EN0 


SUBROUTINE  PLACE  ( IZ , I El , I E2 , IT1 , IT2 , WT ) 


SUBROUTINE  TO  PUT  COEFFICIENTS  FROM  ZONE  EQUATION 
INTO  MASTER  EQUATION 

COMMON  /EON/  EOU( 3,100),I8EGI(10),IEnDI(10), EQUTOT (2,175), 
• SCALE ( 175) .NN.NUNKS.NUNK.NWRT, IUNIT 

NTRMS*IENOl( IZ)-IBEGI(IZ)+1 
I TOT  * IE2-I E 1 + 1 
DO  10  1*1 , ITOT 

I E* !♦ IE  1 -1 
IT.I+IT1-1 
DO  10  d*  1 , NT RMS 
d1*J+lBEGl(IZ)-1 

10  E0UT0T(IT,J1 )* EQUTOT ( IT ,J1 )-fEQU( IE.J)*WT 
RETURN 
END 


SUBROUTINE  POINTF  ( IZ,X1 , Y1 ,X2, Y2,NPT,NBEG.NNEN0, ICTYP.GPAR. 
2 IMAP, INV.X.Y.XPAR.YPAR) 


SUBROUTINE  TO  FIND  POINTS  AT  WHICH  TO  COLLOCATE 

IZ  - ZONE  NUMBER 

XI  + : Y1  - BEGINNING  POINT 

X2  ♦ I Y2  - ENDING  POINT 

NBEG.NNEND  - NUMBER  OF  POINTS  TO  LEAVE  OUT  AT 
BEGINNING, END 

ICTYP  - 1 — STRAIGHT  LINE 

2 — CIRCULAR  ARC 

3 — ELLIPTICAL  ARC 

IMAP  - MAP  COORDINATES  UNLESS  0 

INV  - 1 — MAPPING  from  PARAMETER 

- 2 — MAPPING  TO  parameter 
GPAR ( 1 ) - AXIS  IN  X-DIRECTION( ICTYP.2,3) 

GPAR(2)  - AXIS  IN  Y-DIRECTION( ICTYP*2,3) 

GPAR ( 3 ) - ANGLE  OF  ROTATION  (ICTYP«3) 

GPAR ( 4 ) - XO  ( ICTYP*2 , 3 ) 

GPAR(5)  - YO  ( ICTYP.2,3) 


COMMON  /PI  FAC/  PI .HLFPI , TWOPI 

DIMENSION  X ( 1 ) . Y ( 1 ) , XPAR(1),  YPAR ( 1 ) , GPAR ( 1 ) 
DATA  ICLOCK  /I/ 

GPAR3R=GPAR(3)». 17453292E-1 
xn=npt-i 
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NSTP=NPT-NNEND 
GO  TO  (10,20,30),  ICTYP 
10  FACT1=(X2-X1 )/XN 
FACT2= ( Y2-Y 1 )/XN 
GO  TO  70 
20  XX=X1-GPAR(4) 

YY=Y1-GPAR(0) 

TH1 ■ TANGO ( YY , XX) 

XX=X2-GPAR(4) 

YY  =Y2-GP AR ( 5 ) 

TH2*  TANGO ( Y Y , XX ) 

GO  TO  40 

30  CALL  QUADF  ( GPAR ( 1 ) . GPAR( 2 ) , GPAR3R , GPAR ( 4 ) , GPAR ( 5 ) , XI  ,Y1  , 1 , XX , YY ) 
THI.TANGO(YY.XX) 

CALL  QUADF  (GPAR ( 1 ) ,GPAR(2) , GPAR3R ,GPAR( 4 ) ,GPAR(5) , X2 , Y2 ,1 , XX , YY ) 
TH2=TANG0(YY,XX) 

40  IF  (ICL0CK.E0.-1 ) GO  TO  50 
IF  (TH2.LT.TH1)  TH2*TH2+TW0PI 
GO  TO  60 

50  IF  (TH2.GT.TH1)  TH2=TH2-TW0PI 
60  TH0=.5*(TH1+TH2)+PI 
THO=BRANCH( THO , -PI ) 

TH1  = BRANCH ( TH1  , THO  ) 

TH2  = BRANCH ( TH2 , THO ) 

THETA=(TH2-TH1 )/XN 
70  00  1 10  1 = 1 ,NPT 

IF  (I.LE.NBEG.OR.I.GT.NSTP)  GO  TO  110 

XNN=I-1 

Ilxl-NSEG 

GO  TO  (80,90,90) , ICTYP 
80  X ( 1 1 )=X1+XNN=FACT1 
Y ( 1 1 )=Y1+XNN*FACT2 
GO  TO  110 

90  ANG=TH1 +XNN-THETA 

IF  (ICTYP. EO. 3)  GO  TO  100 
X(  II  )=GPAR(4)+GPAR( 1 ) *COS ( ANG ) 

Y ( 1 1 ) =GP AR ( 5 ) +GPAR ( 1 )=SIN( ANG) 

GO  TO  110 
100  XX=COS(ANG) 

YY  = S I N( ANG) 

CALL  WAP  (XX,YY,CPAR(1 ) ,GPAR(2) , GPAR3R , GPAR ( 4 ) ,GPAR(5), 

* 1 ,X( 1 1 ) , Y( 1 1 ) ) 

110  CONTINUE 

NPT  sNPT-NBEG-NNENO 
IF  (INV.EQ.2)  GO  TO  130 
DO  120  1=1 , NPT 

XPAR( I )=X( I ) 

120  YPAR ( I ) = Y ( I ) 

130  IF  (IMAP.EO.O)  RETURN 
DO  140  1=1 , NPT 

140  CALL  MAPPER  ( XPAR( I ) , YPAR ( I ) , X( I ) , Y( I ) , 12 , INV) 

RETURN 

END 


SUBROUTINE  POWSER(XP,YP.A,B, IZ) 

ROUTINE  TO  COMPUTE  TERMS 
FOR  POWER  SERIES 


COMMON  /P I FAC/  PI.HLFPI.TWOPI 
COMMON  /POWERS/  R( 1 00 ) . TH( 1 00 ) , I SHFT 

COMMON  /SERIES/  ISTYP( 10) , I ZONES , XNOT (10) , YNOT ( 10 ) , NNEG( 10) , 

• NPOS( 10). IZRO(IO) ,IOPRL(10) , 

• IODIM(IO) . IEVRL(10) ,IEVIM(10) 

XX=XP-XNOT( IZ) 

YY= YP-YNOT ( IZ ) 


IPOW=-NNEG( IZ) 


m 


ill 
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NT0T1=NIMEG(  IZ ) +NPOS  ( I Z ) ♦! 

’ 1STRT*1SHFT*1 

R1=5QRTIXX»XX+YY»YY) 

TH1=ATAN2(YY.XX) 

XIPOA»IPOW 

R ( ISTRT)*S0RT(A»A+3»B)»(R1*»IP0W) 

TH( I STRT ) = BRANCH ( AT  AN2 ( B . A ) ♦TH1 ‘XIPOW.-PI) 

DO  10  I *2 ,NT0T 1 

I1«ISHFT+I 

R(  11  )*R( 11-1 )«R1 

10  TH(  1 1 ) * BRANCH ( TH( l 1- 1 )*TH1 ,-PI) 

RETURN 

END 


SUBROUTINE  PRNTMP  (IZ) 

routine  to  print  mapping 

INFORMATION  in  A ZONE 
COMMON  /MAP/  I WAP (5, 2) .PAR AM (5, 2, 5) 

IttPa  1 

PRINT  20.  1Z 

IF  ( IVAPtlZ, 1 ) .CT.O)  GO  TO  10 

PRINT  30 

RETURN 

10  PRINT  10.  IMP,  IMA?(  IZ  . 1 ) . ( PARAM(  IZ.  1 . I ) , 1*1 .5) 

IF  ( IVAP( IZ.2I .LE.C)  RETURN 
I '"•'?«  2 

PRINT  40.  IMP.IKAP1 I Z . 2 ) . ( P ARAM( I Z . 2 , I ) , I » 1 , 5 ) 

RETURN 

20  fuR...-.  1////2BH  MAPPING  FUNCTION  PARAMETERS, 

1 13H  FOR  ZONE  NO. .12) 

30  FORMAT  (23H  .-NO  MAPPING  IN  THIS  ZQNE»») 

40  FORMAT  ( 2 1 H MAPPING  FUNCTION  NC..I2.8H  IN  ZONE/ 

1 13H  MAPPING  TYPE.I2/12H  PARAMETERS SE1 2 . 5/ ) 

END 


C 

C 

C 

C 


SUBROUTINE  PRINTR  (1ZN) 

ROUTINE  TO  PRINT  COEFFICIENTS 
OF  EITHER  A PHI  OR  PSI  FUNCTION 


COMMON  /EON/ 


EQU ( 3 . 1 00 ) , I BEG I ( 1 0 ) , I END 1(10) . EQUTOT (2 , 1 75) , 

• SCALE! 175) .NM.NUNKS.NUNK.NWRT. IUNIT 

COMMON  /SERIES/  ISTYP( 10) . I ZONES , XNOT ( 1 0 ) ,YNOT( 10) ,NNEG(10) . 

• NPOSI 10) . IZRO( 10) , IODRL( 10) , 

. IODI.W(  10)  . IEVRL(  10) , IEVIM(10) 

COMMON  /SOLN/  A( 1 540 1 ) , 1 01 AG( 176) , PH 1(176) 

DATA  ZERO  /O./ 

IZR=IZRO( IZN) 

IOR»ICDRL< IZN) 

IOIaICDlM( IZN) 

I ER* I EVRL ( IZN) 

1 E I * I E V 1 K ( IZN) 

PRINT  40 
NPAaNNEGI IZN) 

NPa-NPA 

NP*RS=NPA+NPOS( IZN)+1 
IE0=2 

IF  (NPA/2.NE. (NPAtI )/2)  IE0=1 

I ND= I SEGI ( IZN) 

00  30  1*1 , NPWRS 

PRTRL=ZERO 

PRTIM.ZERO 


i 
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IF  (NP.EQ.O.AND. IZR.EQ.O)  GO  TO  20 

IF  ((IOR.EO.O.AND.IEO.EQ.1).OR.(I£R.EQ.O.AND.IEO.EQ.2))  GO  TO  10 
PRTRl=PHI ( IND) 

IND.IND+1 

10  IF  ((IOI.EO.O.AND.IEO.E0.1).OR.(IEI.EO.O.ANO.IEO.E0.2))  GO  TO  20 
PRTIM.PHI ( IND) 

InD= INO+1 

20  PRINT  50.  NP.PRTRL.PRTIM 
NP=NP+1 
IE0=3-IE0 
30  CONTINUE 
RETURN 

40  FORMAT  ( 4h  P0W.4X.9HREAL  C0EF.4X.9H I MAG  COEF  ) 

50  FORMAT  (1H  ,I3,2Et3.5) 

END 


SUBROUTINE  PRINTS  (NZ) 

ROUTINE  TO  PRINT  INFORMATION  ABOUT 
EITHER  A PHI  OR  PSI  FUNCTION  IN  A ZONE 

COMMON  /SERIES/  I ST Y P ( 1 0) . I ZONES . XNOT (1 0 ) , YNOT ( 1 0 ) ,NNEG( 1 0 ) . 

* NPOS ( 1 0) , IZR0( 10) , I00RL( 10) . 

• I0DIM( 10) , IEVRL( 10) , IEVIM(10) 

IF  (NZ.GT.5)  GO  TO  10 

PRINT  30,  NZ 
GO  TO  20 
10  NNZsNZ-S 

PRINT  40,  NNZ 

20  n:jn=-nneg(nz) 

IF  (NNEG(NZ) .GT.O)  PRINT  50,  nnn 

IF  (NPCS(NZ).GT.O)  PRINT  60,  NPOS(NZ) 

IF  ( IZRO(NZ) .GT.O)  PRINT  70 

PRINT  BO,  IODRL(NZ) , IODIM(NZ) , IEVRL(NZ) , IEVIM(NZ) 

RETURN 

30  FORMAT  ( //22H  PHI  FUNCTION  FOR  ZONE, 12) 

40  FORMAT  (//22H  PSI  FUNCTION  FOR  ZONE, 12) 

50  FORMAT  (26H  HIGHEST  NEGATIVE  POWER  IS, 14) 

60  FORMAT  (26H  HIGHEST  POSITIVE  POWER  IS, 14) 

70  FORMAT  (51H  ZERO  POWER  INCLUOED(IF  CONSISTENT  WITH  SYMMETRIES)) 
80  FORMAT  ( 34H  COEFFICIENTS  I NC LUDED ( 1 * YES , 0«N0 ) / 

1 1 1 H ODD  POWER, 5X, 10HEVEN  POWER/1 1H  REAL  IMAG.SX, 

2 10HREAL  IMAG/1X, 14, I6.5X, 14, 16) 

END 


{ 


* 

I 
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COMMON  /SERIES/  ISTYP(10),I ZONES , XNOT (10) , YNOT (10) ,NNEG( 10), 

• NPOS(IO), IZRO(IO) ,IODRL(10) . 

* IODIM(  10)  , IEVRL(  10)  , 1C''IM(  10) 

PRINT  20,  I Z 

NZ * I Z 

00  1 0 I * 1 , 2 

IF  ( ISTYP(NZ) . LE . 0 ) GO  TO  10. 

CALL  PRINTS  ( NZ ) 

10  NZ=NZ+5 
RETURN 

20  FORMAT  ( ////32H  SERIES  INFORMATION  FOR  ZONE  NO. ,12) 

ENO 


SUBROUTINE  PRTCOF  (IOPT) 

SUBROUTINE  TO  PRINT  COEFFICIENTS  OF 

PHI  AND  PSI  FUNCTIONS  IN  VARIOUS  ZONES 
IOPT  - INTEGER  CONTAINING  ZONE  NUMBERS 


COMMON  /SERIES/  1ST YP( 1 0) , IZONES , XNOT ( 1 0 ) , YNOT( 10 ) , NNEG( 1 0) 

* NPOS( 10) . IZRO( 10) , IODRL( 10) , 

* 100 IM(10) ■ I E VR  L ( 10) ,1 EV IM( 1 0 ) 

DIMENSION  I ZONE ( 5 ) 

CALL  DECODE  ( IOPT , 5 , I ZONE , NZNZ ) 

IF  (NZNZ.EQ.O)  RETURN 
PRINT  SO 
DO  20  1=1 .NZNZ 

NZ  = I ZONE (NZNZ-I+1 ) 

if  (npos(nz).le.o.and.nneg(nZ).le.O)  go  to  io 

I SER* 1ST YP ( NZ ) 

PRINT  30,  NZ 
PRINT  60.  I SER 
CALL  PRINTR  (NZ) 

10  NZ=NZ+5 

IF  (NPOS(NZ) . LE.O.AND.NNEG(NZ) , LE.O)  GO  TO  20 
ISER*ISTYP(NZ) 

PRINT  AO 
PRINT  60,  ISER 
CALL  PRINTR  (NZ) 

20  CONTINUE 
RETURN 

30  FORMAT  (///9H  ZONE  N0..I2//13H  PHI  FUNCTION  ) 

40  FORMAT  (//13H  PSI  FUNCTION  ) 

50  FORMAT  ( / /// //34H  COEFFICIENTS  IN  PHI, PSI  FUNCTIONS  ) 

60  FORMAT  (12H  SERIES  TYPE, 12/) 

END 


SUBROUTINE  PRTCON  (IOPT) 

ROUTINE  TO  PRINT  EXTRA  PARAMETERS 

IOPT  - VARIABLE  CONTAINING  PARAMETER  NUMBERS 

COMMON  /CONST/  I CON ( 5 ,2) , ISTRT(5) 

COMMON  /SOLN/  A(15401),IDIAG(176).PHI(176) 
DIMENSION  I PRT ( 5 ) 

CALL  DECODE  ( IOPT ,5 , I PRT , NC) 

IF  (NC.EO.O)  GO  TO  40 
00  30  I»1 ,NC 

IC*IPRT(NC-I+1 ) 

IF  ( ICON( IC , 1 ) . EO . 0 . AND. ICQN( IC , 2 ) . EQ • 0 ) GO  TO  30 
CONR  L*0 . 

CONIM.O. 

1 1 - ISTRT ( IC) 

IF  (ICON(IC.I).EO.O)  GO  TO  10 


37 


noooooooonoooooo  oooo 


t 


CONRL=PHI  (II) 

11*11+1 

10  IF  ( IC0N( IC,2) .E0.0)  CO  TO  20 
CON IM*PH 1(11  ) 

20  PRINT  50,  IC , CONRL . CON IM 
30  CONTINUE 
40  RETURN 

50  FORMAT  (///20H  EXTRA  PARAMETER  NO. , 12/1 X , 2E 1 3. 5) 
END 


SUBROUTINE  PRTZON  (IOPT) 

ROUTINE  TO  PRINT  ALL  INFORMATION  IN  ZONE 

IOPT  - VARIABLE  CONTAINING  ZONE  NUMBERS 

COMMON  /SERIES/  ISTYP( 10) , I ZONES , XNOT ( 1 0 ) ,YNOT( 10) ,NNEG( 10) . 

* NP0S( 10) , I ZRO( 10),I0DRL(10), 

* IODIM(IO) ,IEVRL(10),IEVIM(10) 

DIMENSION  I ZONE ( 5 ) 

CALL  CECODE  ( I OPT . 5 , I ZONE . NZNZ ) 

IF  (NZNZ.EO.O)  RETURN 
PRINT  20 
DO  10  1*1 .NZNZ 

IZ=IZCNE(NZNZ-I+1 ) 

IF  (IZ.GT.5.0R.ISTYP(IZ).LE.O.ANO.ISTYP(IZ+5).LE.O)  RETURN 
CALL  PRNTSR  (1Z) 

CALL  PRNTMT  (IZ) 

10  CALL  PRNTMP  ( IZ) 

RETURN 

20  FORMAT  ( 1 HI ) 

END 


SUBROUTINE  PTCOL  ( XP , YP , 1 Z , IOPT ,IREF,WT,IF1,IF2,IT1tIT2) 


SUBROUTINE  TO  COLLOCATE  AT  ONE  POINT  EITHER  FORCES 
OR  DISPLACEMENTS 


IOPT 

IREF 


WT 

IF1.IF2  - 
IT01 , IT02- 


1 —  COLLOCATE  FORCES 

2—  COLLOCATE  DISPLACEMENTS 

1—  NO  REFLECTIONS 

2—  X-AXIS  REFLECTION 

3— -UNIT  CIRCLE  REFLECTION 
WEIGHTING  OF  EQUATIONS 
WHICH  EQUATIONS  TO  SAVE 

WHERE  TO  PLACE  EQUATIONS  TO  BE  SAVED 
IN  MASTER  EQUATIONS! EITHER  1ST  OR  2ND  EQUAT 


COMMON  /MAT/  G(5).ETAM(5) 

C 1 = 1 . 

C 2 = 1 • 

IF  ( IOPT. EQ. 1 ) GO  TO  10 
C 1 * E T AM ( I Z ) 

C2*-1  . 

10  CALL  ZERO  (1,2, 0,0, 0,0) 

CALL  B4PHIP  (XP.YP, IZ, IREF,C2,A,B) 
CALL  PHIPB  (XP.YP.A.B, IZ) 

A*  1 . 

B=0 . 

I F ( IREF.EQ.4)CALL  MAPPR(XP ,YP , A , B , IZ) 
CALL  PHI  (XP.YP, C1*A,-C1*B,IZ) 

GO  TO  (40,20,30,40),  IREF 
20  CALL  PHIB  ( XP , YP , -C2 , 0 . , I Z ) 
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GO  TO  40 

30  CALL  PH1B1  (XP.YP.-C2.0. . IZ) 

40  CALL  PLACE  ( I Z . I F 1 . 1 F2 , I T 1 , 1 T2 . WT ) 

IF  ( IREF.EQ.2.0R. IREF.EQ.3)  GO  TO  50 
CALL  ZERO  (1 ,2. 0,0. 0.0) 

CALL  PHI BR  (XP.VP.C2.0. .12+5) 

CALL  PLACE  (IZ+5,IFt,IF2.IT1,IT2,WT) 
50  RETURN 
END 


SUBROUTINE  PWSERD(XP, YP.A.B, IZ) 


COMPUTER  TERMS  IN  DERIVATIVE 
OF  POWER  SERIES 


COMMON  /PIFAC/  PI .HLFPI.TWOPI 
COMMON  /POWERS/  R( 1 00 ) , TH( 1 00 ) , ISHFT 

COMMON  /SERIES/  ISTYP( 10) , I ZONES , XNOT (10) , YNOT (10) ,NNEG( 10) 

• NPOS(IO), IZRO(IO) .IOORL(IO) . 

• IODIM( 10) , IEVRL(10) ,IEVIM(10) 

NTOT 1 =NNEG ( IZ)+NPOS( IZ)+1 

IPOW*-NNEG( IZ) 

XX*XP“XNOT ( IZ) 

YY*YP-YNOT ( IZ) 

XI POW* I POW 
RHOAB=SQRT(A*A+B*B) 

IF(RHOAB.NE.O. )PHIAB*ATAN2(B,A) 

R1*SQRT(XX*XX+YY*YY) 

TH1 * ATAN2 ( YY , XX ) 

DO  10  1*1, NTOT 1 
11*1 SHFT+ I 
X I POW* I POW 

R( II  )*  RH0AB*XIP0W*(R1**( IPOW-1 ) ) 

TH( 1 1 )* BRANCH ( PH  I A8+ ( XI POW— 1 . )*TH1 ,-PI) 

1 F ( I POW . GE • 0 )G0  TO  10 
R(I1)»-R(I1) 

TH( 1 1 )* BRANCH ( TH( 1 1 )+Pl,-PI) 

10  IP0W*IP0W+1 
RETURN 
END 


SUBROUTINE  OUADF  ( A , B, PHI , XO , YO , XM, VM, IMAP , XQ, VO) 


SUBROUTINE  USED  BY  INVERSE  MAPPINGS  TO  COMPUTE  SQARE  ROOT 
OF  FUNCTION  OF  THE  FORM  (Z«Z  -C*C)  — C IS  REAL 


COMMON  /PIFAC/  PI. HLFPI.TWOPI 
EQUIVALENCE  (PIH.HLFPI),  (PIT.TWOPI) 
DENOMsA+B 
X 1 * XM-XO 
Y 1 * YM-YO 

IF  (DENOM.EQ.O)  GO  TO  50 
DENOM* 1 ./DENOM 
FACT  * 1 • 

IF  (IKAP.EQ.-2)  FACT.-1. 

ICT  = 0 

CS*COS(PHI) 

SN*SIN( PHI ) 

XZ»X1*CS+Y1*SN 

VZ«Y1*CS-X1*SN 

XC«0. 
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XCsSQRT ( AA-BB) 
YC-SQRT(BB-AA) 


VC»0. 

AA>A*A 
BB*B*B 

IF  (AA.GT.BB) 

IF  (BB.GT.AA) 

X1.XZ+XC 
Y1 *YZ+YC 

ri*sqrt(xi*xi«-vi*yi ) 

TH1 >0. 

IF  (XI .EQ.O. .AND.Y1 .E0.0. ) GO  TO  10 
TH1 *ATAN2 ( Y 1 ,X1 ) 

IF  (BB.GT.AA)  TH1 * BRANCH ( TH1 , -3 . *PI ) 

10  X 1 *XZ-XC 
Y1*YZ-YC 

R2*SQRT(X1»XH-Y1*Y1) 

TH2*0 . 

IF  (XI .E0.0. .AND.Y1 .EQ.O. ) GO  TO  20 
TH2* ATAN2 ( Y 1 ,X1 ) 

IF  (BB.GT.AA)  TH2* BRANCH ( TH2 ,-PlH) 

IF  (AA.GT.B8)  TH2* BRANCH ( TH2 ,0 . ) 

20  X2*SQRT(R1 *R2) 

Y2*  . 5*  ( TH1  ♦ TH2 ) 

30  XQ=0ENOM* ( XZ+FACT *X2*C0S( Y2 ) ) 

Y0=0EN0M«(TZ+FACT*X2»SIN( Y2) ) 

XX*XQ»XQ+YQ*YQ 

IF  ( XX . GT . . 9999. OR .IABS(IMAP).E0*2)  RETURN 
IF  (ICT.NE.1)  GO  TO  40 
WRITE  (6.60)  XM,YM 
STOP 

40  FACT »— FACT 
ICT  * 1 
GO  TO  30 

50  XQ  = AT AN2 ( Y 1 ,X1  ) 

IF  (XO.GT.PHI+PI ) XQ.XQ-PIT 
Ir  (XQ.LT. PHI-PI)  XQ.XQ+PIT 
XQ.PHI-XQ 

YQ*AL0G(SQRT(X1*X1-*-Y1*Y1  )/A) 

RETURN 

60  FORMAT  (15H  ERROR  IN  QUADF , 1 OX. 3HXM. , E12 .4, 5X.3H  , E12 .4) 

END 


SUBROUTINE  PWSRD2(XP , YP , A , B. 12) 

ROUTINE  TO  COMPUTE  TERMS  IN 
SECOND  DERIVATIVE  OF  POWER  SERIES 

COMMON  /PI FAC/  PI.HLFPI.TWOPI 
COMMON  /POWERS/  R( 1 00) , TH( 1 00 ) , I SHFT 

COMMON  /SERIES/  ISTYPf 10) , I ZONES, XNOT ( 10) , YNOT( 10) ,NNEG( 10) • 

• NPOS( 10), IZRO(IO) , IODRL( 10). 

* IODIM(10),IEVRL(10),IEV!M(10) 

NT0T1=NNEG(IZ )+NPOS( IZ)+1 

IPOW*-NNEG(IZ) 

XX*XP“XNOT ( IZ) 

Y Y* YP- YNOT ( I Z ) 

R1*SQRT(XX*XX-t-YY*YY) 

TH1aATAN2(YY,XX)  •% 

RHOAB.SQRT(A*A+B*B) 

PHIAB=ATAN2(B.A) 

DO  10  1*1 , NTOT 1 
XIPOW*IPOW-2 
RFACT*IP0W*(IP0W-1 ) 

I1*ISHrT+I 

R( 1 1 )»RH0AB*RFACT*(R1**(IP0W-2) ) 

TH ( 1 1 ) “BRANCH ( PHI AB+X I POW*  TH1 ,-PI) 

10  IPOu/.IPOW+l 


SUBROUTINE  RAOO  (RR , TT , NT RMS) 

ADD  SERIES  OF  COMPLEX  NUMBERS  IN  POLAR  FORM 

COMMON  /POWERS/  R( 1 00 ) . TH( 1 00 ) , I SHFT 
DIMENSION  RR( 1 ) , TT ( 1 ) 

CO  30  1*1 .NT RMS 

I  1 a I SHFT  -f  I 

XX=RR(  1 1 ) »CCS(  TT(H))+R(I1  ) »COS(  TH(  1 1 )) 

YY*RR( 1 1 )«SIN(TT( H ) )+R(1 1 )*SIN(TH( II  ) ) 

IF  (XX.NE.O..OR.YY.NE.O.)  CO  TO  10 
R( II )*0. 

TH( 1 1 )*0. 

GO  TO  30 

10  TH( II )=ATAN2(YY, XX ) 

IF  (ABS(C0S(TH(I1))).LT..1)  GO  TO  20 
R ( 1 1 ) *XX/ COS( TH( II)) 

GO  TO  30 

20  R( II )*YY/SIN(TH( II ) ) 

30  CONTINUE 
RETURN 
END 


SUBROUTINE  RHS  ( I0PT,NPT,X,Y,X1 ,Y1 ,X2 , Y2 . RHS1 0 , RHS20 . SIGO . SIG1 , 

1 RHS1.RHS2) 

ROUTINE  TO  COMPUTE  VALUES  OF  FI  + I F2 
MEANING  OF  IOPT 

1 - CONSTANT  VALUE 

2 - LINEARLY  VARYING  NORMAL  STRESS.  STRAIGHT  BOUNDARY 

3 - LINEARLY  VARYING  SHEAR  STRESS , STRAIGHT  BOUNDARY 

4 - CONSTANT  NORMAL  STRESS,  CIRCULAR , ELLIPTICAL  BOUNDARY 

5 - CONSTANT  SHEAR  STRESS,  CIRCULAR, ELLIPTICAL  BOUNDARY 

COMMON  /PIFAC/  PI.HLFPI.TWOPI 
DIMENSION  X(1),  Y ( 1 ) , RHS1 ( 1 ) , RHS2(1) 

C1*0. 

C2*0 . 

C3*0. 

C4*0 , 

GO  TO  (90,10,20,60,70),  IOPT 
> C1«1  . 

C2»1  . 

GO  TO  30 
I C3*1  . 

C4*  1 . 

) THETAC*ATAN2( Y2-Y1 , X2-X1 ) 

IF  (ABS(THETAC-HLFPI).GT.. 00001)  GO  TO  40 
CSINV«1 . 

GO  TO  90 

) IF  (ABS(THETAC-t-HLFPI).GT.. 00001)  GO  TO  SO 
CSINV*-1 . 

GO  TO  90 

» CSINV-1 ./COS(THETAC) 

GO  TO  90 
I C1«1. 

GO  TO  80 
I C3*1  , 

I CS I NV * 1 . 

I DO  100  1*1 , NPT 

XSODIF* ,5*(X( I )*X( I )-X1*X1 )*SIG1 *CSINV 
VSCDIF* . 5* ( Y( I )*Y(I )-Y1*Y1 )*SIG1*CSINV 
XDIFF-1 X( I )-X1 ).SIG0 
YDIFF*(Y(I)-Y1) *S IGO 

RHS1 ( I ) *RHSlO+C1 *XDI FF+C2»XS0OI F-C3*YDIFF-C4»YS0DI F 
RHS2(  I ) *RHS20+C1  »YDI  FF+C2*YS0DIF+'C3*XDIFF+'C4»XS00I  F 
I CONTINUE 
RETURN 
END 
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SUBROUTINE  ROTEUU  ( IOPT, IZNS, At .A2.A3.A4) 


MULTIPLY  EOUATIONS  IN  GIVEN 
ZONES  BY  COMPLEX  CONSTANTS 


DIMENSION  IZ(5) 

CALL  DECOOE  ( IZNS ,5 . IZ .MAX) 

IE  ( MAX . LE • 0 ) RETURN 
DO  10  1*1 .MAX 

10  CALL  ROTZON  ( IOPT , I Z ( I ) , A 1 , A2 , A3 . A4) 
RETURN 
ENO 


SUBROUTINE  ROTZON(IOPT,IZ.A1 .A2.A3.A4) 


MULTIPLY  SERIES  BY  COMPLEX  NUMBERS 


COMMON  /EON/  EQU( 3, 100). I BEG I (10). 1EN0I ( 10) , EOUTOT (2 , 175) 
• SCALE ( 1 75 ) .NN.NUNKS ,NUNK ,NWRT, I UN IT 

NTRVS=IENDI(  I Z ) - 1 BEGI  ( IZ)  + 1 
DO  30  1*1 , NTRMS 
1 1 *IBEGI ( IZJ  + I-1 
I F ( IOPT . EO • 2 )G0  TO  10 

EOU ( 1 . I ) * EOUTOT (1.11 ) • A1 + EOUTOT (2,11 ) *A2 
10  I F ( IOPT . EQ. 1 )G0  TO  20 

EOU  ( 2 . 1 ) * EOUTOT  ( 1,11)  *A3-f  EOUTOT  (2,11  )*A4 
20  I F ( IOPT .NE.2)E0UTCT ( 1 , II )*0. 

30  IF(IOPT.NE.I) EOUTOT (2,I1)«0. 

11*1 

12*2 

I F ( IOPT .EQ. 1 )I2*1 
IF( IOPT. EQ. 2)11*2 
CALL  PLACE ( IZ, II, 12, 11,12,1.) 

RETURN 

END 


SUBROUTINE  SAVE  (IOPT) 


SUBROUTINE  TO  SAVE  EOUATIONS  STORED  IN  ARRAY  EOUTOT 
IF  IOPT*0  DON'T  DO  ANY  SCALING 


COMMON  /EON/  EOU ( 3 , 1 00 ) , I BEG I ( 1 0 ) , I ENO 1(10), EOUTOT (2,175) 
• SCALE (175) .NN.NUNKS, NUNK.NWRT.IUNIT 

IF  (IOPT.EO.O)  GO  TO  20 
DO  10  1*1,2 

DO  10  J*1 , NUNK 

10  IF  ( ABS ( EQUTOT ( I , J) ) .GT .SCALE(J ) ) SCALE(d) «A8S(EQUT0T( I , J) ) 
20  WRITE  ( I UN  IT)  ( ( EQUTOT ( I , d ) , J*1 ,NN) , I *1 , 2 ) 

NWRT  *NWRT+1 

RETURN 

END 


SUBROUTINE  SCLFIX 


COMMON  /EQN/  E0U(3, 100) , I BEGI ( 10) , I END I ( 10) ,E0UT0T(2, 175) 
* SCA  LE ( 1 75 ) , NN , NUNKS , NUNK , NWRT , I UNIT 

DO  10  1*1 .NUNK 

10  SCALE(I )*1 ./SCALE(I ) 

RETURN 

END 
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subroutine  serinfi iz, isertp, istp, ineg, ipos. izero, icoef 

• IPL.XO.VO.E.XNU) 


SUBROUTINE  TO  SET  VARIABLES  CONTAINING  SERIES 
INFORMATION 

ISERTP  - 1 FOR  PHI  FUNCTION 

2 FOR  PSI  FUNCTION 

ISTP  - EXPANSION  TYPE 

INEG  - NUMBER  OF  NEGATIVE  POWERS 

IPOS  - NUMBER  OF  POSITIVE  POWERS 

IZERO  - 0 IF  ZERO  POWER  EXCLUOED 

IOR  - 0 TO  EXCLUDE  REAL  COEF.  ODD  POWER 

101  - 0 TO  EXCLUDE  IMAG  COEF,  ODD  POWER 

I ER  - 0 TO  EXCLUDE  REAL  COEF,  EVEN  POWER 

IEI  - 0 TO  EXCLUDE  IMAG  COEF,  EVEN  POWER 

I PL  - 1 FOR  PLANE  STRESS 

2 FOR  PLANE  STRAIN 
(XO.YO)  - EXPANSION  POINT 

E.XNU  - ELASTIC  CONSTANTS 

COMMON  /EON/  EQU( 3,100), I BEG I ( • 0) , 1 END I (10), EOUTOT (2,175) 

• SC  A L E ( 1 75 )' , NN . NUNKS . NUNK  , NWR  T , I UN  I T 
COMMON  /MAT/  G(5).£TAM(5) 

COMMON  /SERIES/  1ST YP( 10) , I ZONES ■ XNOT (10), KNOT ( 1 0 ) , NNEG( 1 0 ) , 

• NPOS ( 10) , I ZRO( 10), IODR  L ( 1 0 ) ■ 

• IOOIM( 10) , IEVRL(IO) , IEVIM(IO) 

DIMENSION  I COF ( 4 ) 

EQUIVALENCE  (IC0F(4) , IOR ) , ( ICOF ( 3) , IOI ) , ( ICOF ( 2 ) , I ER ) , 

• (ICOF(t).IEI) 

IZ1.IZ 

IF  ( ISERTP. EQ-2)  IZ1.IZ+5 

izones=maxo( iz, izones) 

ISTYP( IZ1 )»ISTP 
NNEG( IZ 1 ) = I NEG 
NP05(IZ1 )*IPOS 
IZROt IZ 1 ) * IZERO 

CALL  DECODE( IC0EF.4, ICOF.NDUM) 

IOORL( IZ1 )=IOR 
I OD I M ( IZ1 )*IOI 
IEVRL( IZ1 ) = I ER 
I E V IM( 1 Z 1 ) = I El 
XNOT ( I Z 1 )«X0 
YNOT ( I Z 1 )*Y0 
NEVP=I N EG/2+1 PQS/2 
NODP=NEVP 

IF  ( INEG/2.NE. ( INEG+1 )/2)  NOOP»NOOP+1 
IF  (IPOS/2. NE. ( IPOS+t )/2)  N0DP*N00P+1 

IF  ( IZERO. NE.O)  NEVP.NEVP+1 
NZUNK=NEVP*( IER+IEI ) +NODP* ( IOR+IOI ) 

IBEGI ( 121 ) =NUNKS+1 
NUNKS=NUNKS+NZUNK 
IENDI (121 ) =NUNKS 
NUNK  =nunks 

NN*NUNK+1 

G( IZ)=E/(2.*(1 . +XNU ) ) 

GO  TO  ( 10,20) . I PL 
10  ETAM( IZ)=(3.-XNU)/( 1 . +XNU) 

GO  TO  30 

20  ETAM( IZ)=3.-4.*XNU 
30  RETURN 
El 


SUBROUTINE  SERFIL  ( IZ , FI , F2 , A, 8) 


SUBROUTINE  TO  FILL  EQUATIONS  FOR  UNKNOWN  COEFFICIENTS 


COMMON  /EON/  EQU( 3, 100), I BEG I ( 10) , I END I ( TO). E0UT0K2 , 175). 

• SCALE ( 1 7S) , NN.NUNKS , NUNK , NWRT , I UN IT 
COMMON  /POWERS/  R( 100 ) . TH( 100 ) , ISHFT 

COMMON  /SERIES/  ISTYP( 10) . IZONES . XNOT( 10) , VNOT( 10) ,NNEG( 10) . 

• NPOS( 10) . IZRO( 10) , IOORL( 10) • 

• lODIM(IO) , IEVRL(10) , IEVIM( 1 0 ) 

IZR.IZRO(IZ) 

NT  1 *NN£G(  1Z  ) ♦NPOS(  I Z ) ■*•■% 

IOR.IODRL(IZ) 

IOI*IOOIM( IZ) 

UR.IEVRL<  IZ) 

IE!*IEVIM( IZ) 

IPOW— NNEG(IZ) 

IEOal 

1P1«ABS( IPOW) 

IF  (lP1/2.E0.(IP1+l)/2)  I E0*2 

I NO  a 1 

00  30  I«1 ,NT1 

IF  (IZR. EO.O. AND. IPOW. EO.O)  CO  TO  20 
CS«COS(TH(I)) 

SN-SINITH(I)) 

IF  ((IOR.EO.O.ANO.IEO.£Q.1).OR.(IER.EQ.O.AND.1EO.EQ.2))  g0  to  10 
EOU(  1 , 1 NO ) «EOU(  1,IND)*R(I)*(A*CS-B*F1*SN) 

EQ'J  ( 2 , 1 NO ) *EQU  (2.1  NO  ) +R  ( I )»(B»CS+A»F1  • SN ) 

IN0.IN0+1 

10  IF  ((IOI.EO.O.AND.IEO.E0.1).OR.(!EI.EQ.O.AND.IEO.EQ.2))  GO  TO  20 
EQU( 1 , 1ND) «EOU ( 1 ,IND)-F2*R(I)*(B*C5+A»F1»SN) 

EOU( 2 , INO) *EOU ( 2 ,IND)'fF2«R(I)»(A»CS“B*Fl*SN) 

IND.IND+1 
20  IP0W«I?0W*1 
IE0.3-IE0 
30  CONTINUE 
RETURN 
ENO 


subroutine  slvint 


SUBROUTINE  TO  INITIALIZE  VARIABLES  FOR  SOLUTION  ROUTINE 


COMMON  /EQN/  EOU( 3 , 1 00 ) . I BEG I (10), I ENO I ( 10), EQUTOT (2,175) 
• SCALE( 175) .NN.NUNKS, NUNK, NWRT, IUNIT 

COMMON  /SOLN/  A(1 5401), 101 AG (176). PH 1(1 76) 

IOI AG( 1 ) *1 
N * NUNK 

00  10  1*1, NUNK 

IDlAG(If1)*IDIAG(I)+N 
10  N-N-1 

N* IOl AG(NUNK) 

00  20  I>1  , N 

20  A(I)>0. 

00  30  1*1, NUNK 

30  PHI ( 1 ) *0 . 

RETURN 

ENO 


onooo  nonnnonooon 


subroutine  STCOL(XP,YP,N, IZ, IREF, ALPHA, RHS1 ,RHS2,IERRI 


SUBROUTINE  TO  COMPUTE  STRESS  EQUATIONS 
MEANING  OF  IREF 

1 - NO  REFLECTIONS 

2 - X-AXIS  REFLECTION 

3 - UNIT  CIRCLE  REFLECTION 

ALPHA  - ANG..E  OF  NORMAL  TO  SURFACE  ON  WHICH  TO  FIND  STRESS 

DIMENSION  XP( 1 ) , YP( 1 ) ,RHS1 (1 ) .RHS2( 1 ) , ALPHA ( 1 ) 

COMMON  /EQN/  EQU (3, 100) , I8EGI( 10) , I END I (10), EQUTOT (2,175), 

• SCALE ( 175) , NN , NUNKS , NUNK , NWRT , IUNI T 

COMMON  /Pi  FAC/  PI .HLFPI .TWOPI 
DO  20  1=1 ,N 
CALL  ZERO( 0 , 0 , 1 , 2 , 1 , 0) 

ALFA  = Pl *ALPHA( I )/ 1 80  . 

CALL  STRPT(XP( I ) . YP ( I ) . IZ , ALFA. I REF , 1 . , 1 . 2 , 1 , 2 ) 

EQUTOT ( 1 ,NN)sRHS1 ( I ) 

EQUTOT ( 2 , NN ) * RHS2 ( I ) 

I F ( 1 ERR . EQ . 1 ) GO  TO  10 
CALL  SAVE( 1 ) 

GO  TO  20 

10  CALL  EVAL(1 .2.2.XSR. YER) 

PRINT  30,  XER,YER.XP(I).YP(I) 

20  CONTINUE 
RETURN 

30  FORMAT ( 1 H , 2F 1 0 . 4 ,2E 1 0 ,4) 

END 


subroutine  solve 


SUBROUTINE  TO  SOLVE  SIMULTANEOUS  EQUATIONS 


COMMON  /EQN/  EOU( 3 , 1 00 ) , I BEG I ( 1 0) , I END I (10), EQUTOT (2 , 175). 
• SCALE ( 175) , NN, NUNKS , NUNK , NWRT , I UN IT 

COMMON  /SOLN/  A( 1 5401 ),IDIAG( 176), PHI (176) 

EQUIVALENCE  (N.NUNK) 

DIMENSION  T ( 1 76) 

DO  70  1 = 1 , N 

L= IDIAG( I ) 

TEf.lP  = A ( L ) 

LL= I DI AG( !♦! )-1 
1 1 = I“L 

DO  10  U=L,LL 

T( II+J)=A(U) 

10  A ( J ) s A( U ) /TEMP 

PHI(I)=PHI(I)/TEMP 
DO  60  U*1  ,N 

MM=IDIAG( J+1 )-1 
NOFF=LL"MM 
IF  (J-I)  20,60,30 

20  m«io:ag( j)+i-j 

T EMP* A( M ) 

GO  TO  40 
30  M*IDIAS(U) 

TEMP«T( J) 

40  DO  50 

50  A(K)=A(K)-TEMP*A(K+NOFF) 

PHI( J)=PHI (J)-TEMP*PHI (I) 

60  CONTINUE 
70  CONTINUE 
RETURN 
END 
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SUBROUTINE  STITCH  (NZl ,NZ2 , NREF1 . NREF2 . I SOPT , INV, 


c 

1 

INVZN.NPT.XX.YY.FI ,F2,U,V.  IERR) 

c 

c 

SUBROUTINE 

TO  STITCH  TWO  ZONES 

c 

ISOPT 

- 1— STITCH  FORCES  ONLY 

c 

2— STITCH  DISPLACEMENTS  ONLY 

c 

3 — ST  ITCH  FORCES, DISPLACEMENTS 

c 

INV 

- 1 — POINTS  COMING  FROM  PARAMETER  PLANE 

c 

OF  ZONE  INVZN 

c 

c 

IERR 

- 0 IF  NOT  CALCULATING  ERRORS 

c 

COMMON 

/EON/ 

EQU( 3 , 100 ) , IBEGI (10), I ENDI (10), E0UT0T(2  f 175), 

• 

SCALE( 175) .NN.NUNKS.NUNK.NWRT, IUNIT 

COMMON 

/MAT/ 

G ( 5 ) , E T AM  ( 5 ) 

DIMENSION  XX(1), 
IREF1 *NREF1 
IREF2*NREF2 

YY ( 1 ) , FI ( 1 ) , F2( 1 ) , U(1),  V(1) 

I Z 1 *NZ1 
IZ2*NZ2 
IF  (INV. 

EQ.2.0R. 

INVZN. EQ. NZl. AND. INV. EQ.1)  GO  TO  10 

IREF1 *NREF2 
IREF2*NREF1 
IZ1.NZ2 
IZ2.NZ1 
10  I FU1 * 1 
I FU2  = 2 

IF  (ISOPT.E0.1)  IFU2.1 
IF  (ISOPT.EQ.2)  IFU1.2 
00  90  1*1 ,NPT 

GO  TO  (20,30),  INV 
20  XP1 *XX( I ) 

YP1 «YY( I ) 

GO  TO  40 
30  X=XX(I) 

Y* YY( I ) 

40  CALL  MAPPER  ( XP1 , YP1 , X , Y, IZ1 , INV ) 

CALL  MAPPER  ( XP2 , VP2 , X , Y, IZ2 . 2) 

00  90  J*  I FU1  , 1 FU2 
FACT  *— 1 . 

IF(J.EQ.2)FACT*-G(IZ1 )/G( IZ2) 

CALL  ZERO  (0,0,1 ,2.1 .0) 

CALL  PT COL  (XP1.YP1,IZ1,J.IREF1,1.,1,2,1,2) 
IF(IR£F1 .E0.4)CALL  MLTDSP(XP1 ,YP1 ,IZ1 ) 

CALL  PTCOL  (XP2.YP2, IZ2.J, IREF2, FACT, 1,2, 1.2) 
IF(IREF2.EQ.4)CALL  MLTDSP(XP2,YP2,IZ2) 

GO  TO  (50,60),  d 
50  EOUTOT ( 1 , NN ) * F 1 ( I ) 

E0UT0T(2.NN)*F2( I) 

GO  TO  70 

60  EOUTOT ( 1 ,NN)*U( I ) 

EOUTOT ( 2 , NN ) * V ( I ) 

70  IF  ( I ERR . EO . 0 ) GO  TO  80 
CALL  EVAL  ( 1 , 2 , 2 > XER , YER ) 

PRINT  100,  XER.YER.XP1 ,YP1 .XP2.YP2 
GO  TO  90 

80  CALL  SAVE  (1  ) 

90  CONTINUE 
RETURN 

100  FORMAT  ( 1 H , 2F10.4.4E10.4) 

END 


SUBROUTINE  STRPT(XP.YP,IZ,ALFA,IREF,WT,IF1 .IF2.IT1 ,IT2) 

c 

c 

C SUBROUTINE  TO  COLLOCATE  STRESSES  AT  A POINT 

C XP  ♦ I VP  - POINT  AT  WHICH  TO  COLLOCATE 
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IZ  - ZONE  NUMBER 

ALFA  - ANGLE  OF  NORMAL  TO  SURFACE  AT  WHICH  STRESS  I 

IREF  - REFLECTION  OPTION 

1 — NO  REFLECTIONS 

2 — X-AXIS  REFLECTION 

3 — UNIT  CIRCLE  REFLECTION 

WT  - WEIGHTING  OF  EQUATIONS 


ALPHA. -ALFA 
CSA.COS( 2 . .ALPHA) 

SNA.S IN ( 2 . * ALPHA ) 

CALL  ZERO( 1,2. 0,0, 0.0) 

I F ( I FI .GT . 1 )G0  TO  10 
CALL  MAPPR(XP, YP.A.B, IZ) 

DEN0M.1 . / ( A»A+B*B) 

AA.2. *A» DENOM 
BB=-2. *B» DENOM 

call  phipb(xp,yp. a a , -bb , I z ) 

CALL  ZERO( 2 , 2 , 0 , 0 . 0 , 0 ) 

10  CALL  B4PH0P(XP,YP,IZ,CSA,SNA,IREF,A.B) 
CALL  PHIPPBfXP. YP,A,-B, IZ) 

CALL  B4PHP(XP, YP.CSA.SNA, IZ, IREF ,A,B) 
CALL  PHIPB(XP,YP.A.-B. IZ) 

CALL  MAP PR (XP,YP,XPR,YPR,IZ) 

CSA.C0S(2. .ALPHA) 

SNA. SIN(2. .ALPHA) 

DENOM. 1 ./(XPR.XPR+YPR.YPR) 

A« DENOM* (CSA*XPR+SNA»YPR) 

B*  DENOM  * ( SNA*  XPR-CS A*  YPR ) 

GO  TO  (40,20,30) , IREF 
20  CALL  PHI8P1(XP,YP,A,-B,IZ) 

GO  TO  40 

30  XT.XP*XP-YP*YP 
YT=2.*XP*YP 
DENOM: 1 ./(XT*XT+YT*YT) 

AA.-DENOM*  (A»XT+-B*YT) 

BB*  DENOM* (A*YT-B»XT) 

CALL  PHIPiB(XP,YP,AA,-BB.IZ) 

40  CALL  PLACE(IZ.IF1.IF2,IT1,IT2,WT) 

I F ( IREF.NE. 1 )G0  TO  50 
A*-A 

gr_0 

CALL  Z£RO( 1 ,2, 0,0, 0,0) 

CALL  PHIPB(XP,YP,A.-B, IZ+5) 

call  PLACE(IZ+5,IF1 , IF2.IT I .IT2.WT) 

50  RETURN 
END 


SUBROUTINE  STRRES  (IZ.NPT, IREF, ALFA, X.Y.XP.YP) 


subroutine  to  compute  and  print  stresses  from  solution 

IREF  - 1 FOR  NO  REFLECTIONS 

2 FOR  X-AXIS  REFLECTION 

3 FOR  UNIT  CIRCLE  REFLECTION 

ALFA  - ANGLES  OF  NORMAL  TO  SURFACES  ON  WHICH  ST 

C ARE  CALCULATED 


DIMENSION  X(1),  V(1),  XP(1),  YP(1),  ALFA(1) 

COMMON  /PIFAC/  PI , HLFPI , TWOP I 

FACT.1B0./PI 

PRINT  20 

1,0  10  I«1  , NPT 

CALL  ZERO  (0.0,1 ,2,3.0) 

CALL  STRPT  <XP(I). YP(I),IZ,ALFA(I), IREF, 1., 1,2, 1,2) 


nnnnn  non 


CALL  EVAL  (1,2,1 .SIG.TAU) 

CALL  ZERO(0,0. 1 ,2,3.0) 

CALL  STRPT(XP(I),YP(I),IZ.ALFA(I)*HLFPI,IREF,1.,1.1,1.1) 
CALL  EVAL( 1 . 1 . 1 .SIG90.SOUM) 

ALF»ALFA( I )»FACT 

10  PRINT  30,  SIG.SIG90,TAU,ALF.X(I),Y(I),XP(I),YP(I) 

RETURN 

20  FORMAT  (////23HOEVALUATION  OF  STRESSES  / 7X , 7HSTRESS1 . 

* 6X.7HSTRESS2, IX, 12HSHEAR  STRESS , SX , 5HANGLE , 

• 9X , 1 HX , 1 2X , 1 HY , 8X , 8HX ( PARAM) , 5X , 8HY ( PAR AM ) ) 

30  FORMAT  (1X,3E13.5,F10.3,4E13.5) 

END 


FUNCTION  TANGO  (XI. X2) 
TANGO=BRANCH( ATAN2(X1 ,X2) ,0. ) 
RETURN 
END 


SUBROUTINE  UNIT(I) 

SUBROUTINE  TO  SET  UP  NEW  TAPE  FOR  WRITING  EQUATIONS 

COMMON  /EQN/  EQU( 3,100),IBEGI(10),IENDI(10), EQUTOT (2,175), 
• SCALE ( 175) , NN , NUNKS , NUNK , NWRT , I UN  I T 

IUNIT  = 1 
NWRT  *0 
RETURN 
END 


SUBROUTINE  WRTDSP(NF3,NFD, IZ.IREF) 

ROUTINE  TO  WRITE  DISPLACEMENTS  AT  INTERFACE 
OF  MMC  AND  FINITE  ELEMENT  REPRESENTATIONS  FOR 
BACK SUB ST  I TUT  ION  INTO  FINITE  ELEMENT  EQUATIONS 

COMMON  /MAT/  G(5) ,ETAM(5) 

COMMON  /MMCFE/  X ( 1 50 ) , Y ( 1 50 ) , FORHS ( 300 ) , EQNS( 1 20 , 1 76 ) , 
* KUV( 120.2) 

DIMENSION  FD( 1 ) 

EQUIVALENCE  ( FD( 1 ) , FORHS( 1 ) ) 

REWIND  NFB 
REWIND  NFD 

READ  (NFB)  N1.N2.ID0FST 
NPT=N2/2 

READ  (NFB)  ( X ( I ) , Y ( I ) , I « 1 , NPT ) 

A»1  . 

B = 0. 

DO  100  1=1 , NPT 

CALL  ZERO(0,0, 1 ,2,1 .0) 

CALL  MAPPER(XP, YP,X( I ),Y(I),IZ,2) 

CALL  PTC0L(XP,YP,IZ,2. IREF.1 . ,1 ,2,1 ,2) 

IF(IREF.EQ.A)CALL  M APPR ( XP , YP , A , B , I Z ) 

FACT  = 1 ./( ( A*A  + B*B)*2. *G( IZ) ) 

CALL  EVAL ( 1 , 2 , 1 , V 1 , V2 ) 

FD(2*I-1 ) a ( VI  * A-V2* 8 ) * F ACT 
100  FD( I »2  ) * ( VI •B+V2*A)*FACT 

WRITE  (NFD)  ( FD( I ) , I = 1 ,N2) 

END  FILE  NFD 

RETURN 

END 
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SUBROUTINE  ZETN(XP,YP, IZ) 


ROUTINE  TO  CALCULATE  TERMS 
IN  EXPANSION  OF  PHI.  PSI 

COMMON  /SERIES/  ISTYP( 10) . I ZONES . XNOT ( 1 0 ) , YNOT( 10) ,NNEG( 10) . 

• NPOS(lO), IZRO(IO) ,IODRL(10) . 

• IOOIM(TO) . IEVRL( 10) . IEVIM(IO) 

ITYP=ISTYP( IZ) 

CO  TO  ( 10,20) . ITYP 
10  Ail. 

B=  0 . 

GO  TO  40 
20  YY=YP+1 . 

DEN= 1 ./<XP-XP+YY«YY) 

As  XP-OEN 
B = -YYOEN 

AO  CALL  PO'ftSER ( XP.YP.A.B.IZ) 

RETURN 

END 


SUBROUTINE  ZERO  ( I E 1 . I E2 . I SI , I S2 , IOPT , 1 1 NPT ) 


SUBROUTINE  TO  ZERO  OUT  ZONE  AND  MASTER  EQUATIONS 


IE1.IE2  -BEGINNING  AND  ENDING  EQUATION  TO  ZERO  OUT 
(ZONE  EQUATIONS) 

IS1.IS2  -BEGINNING  and  ending  equation  to  ZERO  OUT 
(MASTER  EOUATIONS) 

IOPT  -OPTION  FOR  ZEROING  OUT  MASTER  EQUATIONS 

1—  ZERO  UP  TO  CONSTANTS 

2—  ZERO  TO  TOTAL  UNKNOWNS 

3—  ZERO  TO  TOTAL  UNKNOWNS  ( RHS) 

4“ ZERO  TO  UNKNOWN  IINPT 


COMMON  /EON/  EQU (3,100) .1  BEG  I (lO).IENDI(IO), EQUTOT (2,175). 
• SC A LE ( 175) .NN.NUNKS.NUNK.NWRT.IUNIT 

IF  ( IE2-IE1 .LT.O.CR. ( I E2. EQ. 0 . AND. IE1 .EQ.O) ) GO  TO  20 
DO  10  I = I E 1 . IE2 
DO  1 0 J*1 , 100 

10  EQU ( I . J)=0. 

20  IF  ( IS2-IS1 .LT.O. OR. (IS2. EQ.O. AND. IS1 .EQ.O) ) RETURN 
NT.NUNKS 

IF  (I0PT.EQ.2)  NT*NUNK 
IF  (ICPT.EQ.3)  NT*NN 
IF  (I0PT.EQ.4)  NT*IINPT 
DO  30  I » I S 1 . IS2 
DO  3C  J= 1 .NT 
30  EOUTOT ( I ,d)*0. 

RETURN 

END 


SUBROUTINE  ZETN1 (XP. YP, IZ) 

ROUTINE  TO  COMPUTE  TERMS  IN  FIRST 
DERIVATIVE  OF  EXPANSION  FOR  PHI,  PSI 

DIMENSION  RT( 100) ,TT( 100) 

COMMON  /POWERS/  R ( 1 00 ) , TH( 1 00 ) . I SHFT 

COMMON  /SERIES/  I ST YP ( 1 0 ) . I ZONES , XNOT ( 1 0 ) , YNOT ( 10) ,NNEG( 10 ) . 
> NPOS( 10) . IZRO( 10) , IOORL( 10) . 

. IODIM( 10) , IEVRL( 10) . IEVIM(IO) 

NT0T1 »NNEG( IZ)+NPOS( I Z )*1 
ITYP«ISTYP( IZ) 


4 '•»*  • in*,.*  *•  r«-- 


GO  TO  (10,20). ITYP 
10  A* 1 , 

6*0. 

GO  TO  40 
20  YP1 *YP*1 . 

DAB=1  ./(XP=XP+YP1*YP1  ) 
DCDzDAQ*  DAB 
A*  XP*DA3 
B=-YP1 » DAB 

C=0CD*(YP1*YP1-XP-XP) 
D=DC0*XP» YP1 *2 . 

40  CAUL  PNSERD'XP.yP.A.B, IZ) 
I F ( ITYP.EQ. 1 )G0  TO  60 
DO  50  1=1 , NTOT 1 
RT(I)»R(I) 

50  TT(I)»TH(I) 

CALL  POWSERtXP, YP.C.D,  IZ) 
CALL  RADD(RT,TT,NT0T1 ) 

60  RETURN 
END 


SUBROUTINE  ZETN2(XP,YP,IZ) 


routine  to  compute  terms  IN  SECOND 
DERIVATIVE  OF  EXPANSION  FOR  PHI,  PSI 


DIMENSION  RT( 100)  ,TT( 100) 

COMMON  /PI  FAC/  PI.HLFPI.TWOPI 
COMMON  /POWERS/  R( 100) . TH( 100) , ISHFT 

COMMON  /SERIES/  I ST Y P ( 1 0 ) , I ZONES , XNOT ( 1 0 ) , YNOT ( 1 0 ) , NNEG ( 1 0 ) 

• NPOS( 10) , IZRO( 10) , IOORL( 10) , 

* IODIM( 10) , IEVRL( 10) , IEVIM( 10) 

NT0T1=NNEG( IZ ) +NPOS ( I Z ) +1 

ITYP=ISTYP( IZ) 

GO  TO  ( 10,20) , ITYP 
10  A* 1 . 

B=0. 

GO  TO  30 
20  YP1 *YP+1 . 

DEN=1 ,/(XP*XP+YP1*YP1 ) 

As  XP-DEN 
Bs-YP1 »DEN 

30  CALL  PWSRD2 ( XP , YP , A , B , IZ ) 

IF( ITYP.EQ. 1 )G0  TO  60 
DO  40  1=1 , NTOT 1 
RT(I).R(I) 

40  TT(I).TH(I) 

CALL  ZETN1 (XP, YP, IZ) 

DENs2.*S0RT(DEN) 

TTHs-ATAN2(YP1 ,XP)+PI 
DO  50  1=1 , NTOT 1 
R ( I ) sR( I ) » DEN 

50  TH(I)=BRANCH(TH(I)+TTH,-PI) 

CALL  RADD(RT,TT,NT0T1 ) 

60  RETURN 
END 
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